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CHAPTER  1.  INTRODUCTION 


The  imposition  of  a  sinusoidal  oscillation,  of  amplitude  a^  and 
angular  frequency  co,  on  the  mean  flow  through  a  circular  pipe  causes  a 
periodically  varying  pressure  and  velocity  field.  If  the  amplitude  is 
small  enough  a  linear  response  is  obtained  for  which 


=  dp  + 

dx  dx  o 

u  =  u(y)  +  a 


dp_ 

dx 


cos 


(wt)  , 

cos  (tat  +  0  (y) ) , 


(1.1) 

(1.2) 


with  x  being  the  distance  in  the  flow  direction,  y,  the  distance  from 
the  wall  and  t,  the  time.  The  time  mean  pressure  gradient  and  the  time 
mean  velocity  of  the  undisturbed  flow  are  designated  by  and  u(y);  and 


the  amplitudes  of  the  oscillations,  by  a 


d£ 

dx 


and  a 


u(y) 


These 


oscillations  in  the  velocity  field  give  rise  to  oscillations  in  the 
shear  stress  at  the  wall, 


T 

W 


t  +  a  t 
w  o  w 


cos  (wt  +0  ) , 

w 


(1.3) 


with  0  =  lim  0  (y) .  This  thesis  presents  the  results  of  experiments 

w  y-*0 


with  a  fully  developed  pipe  flow  in  which 


w 


and  0  are  related  to 

T 


dx 


w 


The  experiments  were  conducted  at  high  enough  frequencies  so  that 
over  most  of  the  flow  field  inertia  is  dominant;  that  is, 


P 


9u  _  3p 
9 1  9x 


(1.4) 


A  solution  of  (1.4)  for  a  sinusoidally  varying  pressure  gradient  gives 


(1.5) 


2 


up 


A 

u 

d£ 

dx 

and 


0  = 
u 


(1.6) 


Equation  (1.5)  reveals  that,  for  this  high  frequncy  range,  even  small 
oscillations  in  the  volumetric  flow  can  cause  oscillations  in  the 
pressure  gradient  many  times  larger  than  the  mean  pressure  gradient. 

Close  to  the  wall  turbulent  and  viscous  stresses  become  important 
so  that  u  is  related  to  the  pressure  gradient  through  the  eqation 


3u  3p  ,  a2u  ,  3t^' 

P  3t  "  "  3y2  ** 


(1.7) 


The  principal  theoretical  problem  in  predicting  how  the  oscillations 

in  the  wall  shear  stress  are  related  to  the  oscillations  in  the 

pressure  gradient  is  the  determination  of  how  an  oscillating  pressure 

gradient  affects  the  turbulence.  The  experiments  therefore  provide  a 

test  of  closure  methods  that  predict  x^  close  to  a  wall. 

The  motivation  for  this  work  has  come  from  studies  of  the  wall  shear 

stress  variation  along  a  solid  wavy  surface  over  which  a  turbulent  fluid 

is  flowing.  It  was  found,  in  these  experiments,  that  the  amplitude  and 

phase  angle  characterizing  the  periodic  variation  of  the  shear  stress 

along  the  wavy  surface  depend  on  the  wave  number,  a  =  2t t/A,  and  that 

this  dependency  could  be  scaled  by  using  the  kinematic  viscosity,  v, 

and  the  friction  velocity,  u*.  Of  most  interest  was  the  discovery  of 

a  sharp  change  in  phase  angle  with  increasing  dimensionless  wave  number 

+ 

at  a 


.0005 


.001. 


The  interpretation  of  these  results  depends  on  the  prediction  of 
the  wave-induced  variation  of  the  Reynolds  stress  in  the  viscous  wall 
region  (0  ^  y+  ^  40) .  The  turbulence  modeling  differs  from  that  for  flow 
over  a  flat  plate  in  that  the  waves  induce  a  periodic  variation  of  the 
streamline  curvature  and  of  the  pressure  gradient  at  the  surface. 
Thorsness,  Morrisroe  and  Hanratty  (1978)  and  Abrams  and  Hanratty  (1984) 
have  suggested  that  the  observed  behavior  of  the  phase  angle  is  primarily 
associated  with  the  influence  of  wave-induced  pressure  gradient.  They 
argue  that  at  small  a+  an  equilibrium  situation  exists  whereby  the  flow 
behaves  similarly  to  that  observed  in  gradually  expanding  or  converging 
channels.  Regions  of  the  wave  surface  with  a  favorable  pressure  gradient 
experience  a  damping  of  the  turbulence  in  the  viscous  wall  region;  regions 
with  an  unfavorable  pressure  gradient  experience  an  enhancement  of  the 
turbulence.  As  a+  increases  the  pressure  gradient  varies  so  rapidly 
along  the  wave  surface  that  the  turbulence  in  the  viscous  wall  region 
does  not  respond  immediately.  There  is  a  relaxation  from  an  equilibrium 
turbulence  to  a  frozen  turbulence,  manifested  by  the  sharp  change  in  the 
phase  angle  with  increasing  a+. 

The  principal  goal  of  the  present  experiments  was  to  investigate 

the  above  interpretation  by  carrying  out  studies  of  the  influence  of  a 

rapidly  varying  pressure  gradient  on  the  behavior  of  the  viscous  wall 

region,  without  having  to  deal  with  additional  complications  associated 

with  streamline  curvature.  Of  particular  interest  is  the  determination 

of  whether  the  variation  of  0  with  shows  the  same  relaxation 

x 

w 

phenomenon  observed  for  flow  over  wavy  surfaces. 

In  order  to  meet  the  above  objective  it  was  necessary  to  design  the 
experiment  so  that  was  of  the  same  order  as  the  median  frequency  of 


the  turbulent  velocity  fluctuations  in  the  vicinity  of  the  wall 
(~  .01  x  2it).  At  the  same  time,  it  was  desirable  to  avoid  the  use  of 
imposed  oscillations  of  unrealistically  high  frequency.  This  was  done 
by  using  a  system  which  has  low  frequency  turbulence;  i.e.,  water  flow 
in  a  19.4  cm  pipe.  Sinusoidal  oscillations  with  frequencies  of  0.325  Hz 
and  0.625  Hz  were  introduced  into  the  system  with  a  plunger  type  pump. 
The  amplitudes  of  the  oscillations  were  made  small  enough  that  the 
pressure  gradient  and  wall  shear  stress  varied  sinusoidally  with  time. 
The  shear  stress  variation  at  the  wall  was  measured  with  electrochemical 
probes  mounted  flush  with  the  wall.  It  was  not  possible  to  design  this 
experiment  so  that  the  frequency  response  of  the  probe  did  not  have  to 
be  taken  into  account.  An  advantage  of  using  the  electrochemical 
method  is  that  an  analytical ,  rather  than  an  experimental,  scheme  can  be 
used  to  correct  for  frequency  response. 

A  number  of  previous  investigators  have  studied  the  influence  of 
controlled  flow  oscillations  on  turbulent  flows.  An  excellent  review 
has  recently  been  presented  by  Carr  (1981) .  The  work  presented  in  this 
thesis  differs  in  that  attention  is  focused  on  the  wall  shear  stress 
variation  and  on  high  values  of  io+.  Previous  investigators  have  found 
that  the  influence  of  the  oscillations  on  the  mean  velocity  profile  is 
small  and  that  their  effect  on  turbulence  is  felt  principally  in  the 
region  close  to  the  wall.  For  this  reason,  there  exists  a  need  for 
measurements  of  the  wall  shear  stress  oscillations.  In  most  previous 
studies  accurate  velocity  measurements  were  not  made  close  enough  to 
the  wall  to  determine  the  time  varying  velocity  gradient  at  the  wall. 

Recently,  Tu  and  Ramaprian  (1983)  reported  on  direct  measurements  of 
the  wall  shear  stress  with  flush  mounted  wall  heat  transfer  probes. 
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These  were  at  values  of  of*”  too  low  to  observe  the  relaxation  phenomenon 
in  which  we  are  interested.  In  addition,  there  are  uncertainties 
about  these  results  since  no  attempt  was  made  to  take  into  account  the 
frequency  response  of  the  probe. 

The  measurements  of  the  velocity  gradient  at  the  wall  in  this  study 
confirms  previous  findings  by  a  number  of  researchers  that  the  mean  flow 
field  is  actually  not  affected  by  the  imposed  flow  oscillation.  The 
measurements  in  this  work  have  also  verified  the  assumption  that  there 
exists  a  linear  response  of  the  wall  shear  stress  to  the  imposition  of 
a  small  amplitude  flow  oscillation.  Moreover,  it  is  found  that  the 
relative  phase  and  amplitude  which  characterize  this  response  can  be 
correlated  by  a  single  parameter  m‘,  which  combines  the  effects  of 
frequency  and  Reynolds  number  in  the  range  of  variables  studied. 

The  most  striking  feature  of  the  results  is  a  relaxation  effect, 
similar  to  what  has  been  observed  for  flow  over  a  wavy  surface,  whereby 
the  phase  angle  characterizing  the  temporal  variation  of  the  wall  shear 
stress  undergoes  a  sharp  change  over  a  rather  narrow  range  of  At 

(/  larger  than  the  median  frequency  of  the  turbulence,  there  appears 
to  be  an  interaction  between  the  imposed  flow  oscillation  and  turbulence 
fluctuations  in  the  viscous  sublayer,  which  is  not  described  by  present 


theories  of  turbulence. 
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CHAPTER  2.  LITERATURE  SURVEY 


Unsteady  laminar  flows  in  simple  cases  were  solved  analytically 
or  asymptotically.  These  solutions  are  summarized  by  Schlichting  (1979). 
G.  Stokes  (1851)  and  later  Lord  Rayleigh  (1911)  considered  flow  over 
an  infinite  flat  plate  which  executes  linear  harmonic  oscillations 
parallel  to  itself — Stokes  second  problem.  Th.Sexl  (1930)  and  S.  Uchida 
(1956)  treated  the  case  of  flow  through  a  pipe  under  the  influence  of  a 


found 


periodic  pressure  gradient.  An  unique 


to  characterize  periodic  laminar  flow  in  the  pipe,  where  w  is  angular 
frequency,  v,  kinematic  viscosity  and  r  ,  radius  of  the  pipe.  At  very 
low  values  of  ft,  the  velocity  profile  of  the  oscillation  component  is 
in  phase  with  the  imposed  variation  of  pressure  gradient  and  behaves 
like  that  in  steady  state  flow  with  the  same  mean  flow  rate.  At  very 
large  values  of  ft,  the  oscillation  flow  field  moves  like  a  solid  body 
in  the  core  region  of  the  pipe  and  varies  rapidly  close  to  the  wall. 

For  this  case  the  phase  of  the  oscillation  velocity  lags  of  the  pressure 
gradient  by  90°  in  the  center  of  the  pipe,  but  only  by  45°  at  the  neigh¬ 
borhood  of  the  wall.  This  indicates  that  for  high  frequencies  viscous 
effect  plays  an  important  role  in  the  wall  region  and  that  the  outer  flow 
can  be  considered  as  inviscid. 

Because  the  complexities  and  difficulties  encountered  in  experi¬ 
mental  studies  of  unsteady  turbulent  flow,  researches  in  this  area 
were  rarely  seen  in  literature  a  decade  ago.  Interest  in  this  area 
has  been  growing  rapidly  in  recent  years.  Partially,  this  is  because 
results  from  such  studies  are  needed  in  solving  many  urgent  engineering 
problems.  More  importantly, this  recent  interest  can  be  associated  with 
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the  development  of  instrumentation  and  the  availability  of  new  techniques 
of  data  acquisition  and  processing.  Researchers  have  been  mainly  con¬ 
centrated  on  three  subjects:  Oscillating  flow  at  transitional  Reynolds 
number,  turbulent  boundary  flow  over  flat  plate  under  periodic  exitation 
of  free  stream  flow,  and  turbulent  oscillating  pipe  flow. 

In  this  chapter,  a  brief  survey  of  studies  on  these  subjects  will 
first  be  presented.  Then,  studies  in  a  relevant  field,  turbulent  flow 
over  wavy  surfaces,  will  be  discussed.  Finally,  what  has  been  learned 
from  these  studies,  and  what  are  still  unclear  or  controversial  points 
will  be  summarized. 

2.1  Oscillating  Flow  at  Transitional  Reynolds  Numbers 

Experience  with  steady  flows  subjected  to  spatial  pressure  gradients 
has  shown  that  pressure  gradients  can  have  a  significant  effect  on  the 
critical  Reynolds  number.  In  an  unsteady  flow,  one  would  intuitively 
expect  that  the  effects  of  unsteadyness  on  the  flow  structure  will  be 
stronger  on  flows  in  the  neighborhood  of  transition  than  on  flows  at 
very  large  Reynolds  numbers.  This  is  because  the  transition  process  can 
be  sensitive  to  the  strong  acceleration/deceleration  occurring  in  the 
unsteady  flow. 

Gerrard  (1971)  studied  pulsating  flow  at  a  mean  flow  Reynolds  number 
of  3370  in  a  pipe  of  diameter  3.81  cm  with  an  amplitude  of  oscillation 
0.57  of  mean  and  a  period  of  12.2  s.  A  sheet  of  dye  was  produced  by 
electrolysis  from  a  fine  wire  stretched  across  the  diameter.  Cine- 
photography  was  used  to  record  the  motion  of  dyes  and  to  find  the  velocity 
profile.  The  most  significant  information  obtained  from  the  experiment 
was  that  the  flow  was  laminarized  during  the  acceleration  phase  of  flow 
cycle. 
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Ramaprian  &  Tu  (1980)  found  a  similar  phenomenon  in  their  experiments, 
in  which  oil  flow  was  used  in  a  circular  pipe  at  a  mean  Reynolds  number 
of  about  2100.  Two  frequencies  (1.75  Hz  and  0.057  Hz)  were  studied  with 
an  amplitude  30%  of  the  mean  flow.  Instantaneous  velocities  were  measured 
by  using  laser-Doppler  anemometry  (LDA) .  They  observed  that  flow 
oscillation  increased  the  critical  Reynolds  number  and,  under  certain 
conditions,  even  brought  about  laminarization  of  the  flow. 

Shemer  &  Wygnaski  (1981)  approached  the  problem  differently.  They 
compared  the  responses  of  laminar  and  turbulent  air  flow  to  an  imposed 
oscillation  at  the  same  mean  Reynolds  number  4,000.  The  undisturbed 
flow  was  laminar,  and  turbulence  was  triggered  artificially  by  inserting 
a  cylinder  into  the  pipe.  The  hot  wire  anemometer  was  used  to  measure 
the  instantaneous  velocity.  They  found  that  turbulent  flow  responded 
to  an  imposed  pressure  variation  differently  from  laminar  flow  for  the 
same  Reynolds  number  and  the  same  pressure  variation.  They  argued  that 
the  frequency  parameter  no  longer  controls  the  flow  because  the 
relevant  viscosity  is  no  longer  v  but  turbulent  eddy  viscosity  v  , 
which  is  orders  of  magnitude  larger  than  v.  The  effective  Q  is  thus 
much  lower.  They  also  observed  a  partial  laminarization  of  the 
trubulent  flow  during  the  acceleration  of  the  flow.  The  time  averaged 
Reynolds  stress  -ufvf  and  the  phase  averaged  Reynolds  stress  - < uTvT  > 
were  measured  with  an  x-wire.  It  was  observed  that  -uTv!  was  indepen¬ 
dent  of  the  flow  oscillation  and  that  the  phase  of  -<uTvf>  lagged 
behind  the  pressure  variation. 

There  have  been  some  efforts  to  consider  the  stability  of  time 
periodic  flows.  Davis  (.1976)  reviewed  research  works  in  this  area. 

The  term  unstable  for  a  basic  periodic  state  is  defined  by  Davis  as 
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a  situation  for  which  a  disturbance  experiences  net  growth  over  each 
modulation  cycle.  A  state  for  which  every  disturbance  decays  at  every 
instant  is  defined  as  stable.  The  most  interesting  result  coming  from 
these  studies  is  that  there  may  be  a  state  which  is  neither  unstable 
nor  stable;  i.e.,  the  basic  state  is  subject  to  a  disturbance  that  grows 
during  part  of  the  cycle,  that  attains  an  appreciable  amplitude,  that 
decays  during  another  part  of  the  cycle.  This  type  flow  is  called 
transiently  stable  by  Davis.  Both  theoretical  and  experimental  works 
indicate  that  Stokes  layers  are  stable  to  small  disturbances  but  are 
apparently  unstable  to  disturbances  of  large  amplitude.  For  a  fully 
developed  pipe  flow,  Sarphaya  (1966)  found  that  the  critical  Reynolds 
number  increases  with  increasing  amplitude  until  it  reaches  a  peak  value, 
after  which  it  falls  to  values  below  the  critical  Reynodls  number  for 
steady  flow.  More  works  need  to  be  done  in  this  area. 

2.2  Periodic  Turbulent  Boundary  Flow  Over  Flat  Plate 

The  first  experimental  study  on  periodic  turbulent  boundary  flow 
over  a  flat  plate  was  done  by  Karlsson  in  1959.  Amplitudes  of  the 
sinusoidal  oscillation  up  to  34%  of  the  mean  velocity  and  frequencies 
ranging  from  0  to  48  cycles/sec  were  used  in  a  boundary  layer  wind 
tunnel  with  zero  pressure  gradient.  The  hot  wire  was  used  to  measure 
the  mean  velocity,  the  amplitude  of  in- and  out-of  phase  components 
of  the  first  harmonic  of  the  periodic  oscillation,  and  the  intensity 
of  higher  harmonics  and  turbulence.  One  of  the  most  remarkable 
results  of  his  work  is  that  effects  due  to  non-linear  interactions 
between  turbulence  and  the  oscillation,  even  for  oscillation  amplitude 
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as  large  as  34%  of  the  free  stream  velocity,  are  so  small  that  the 
effect  of  the  oscillation  on  the  mean  flow  is  negligible. 

The  group  at  ONERA/CERT  in  France  has  done  extensive  studies  on 
a  turbulent  boundary  layer  subjected  to  sinusoidal  free-stream 
oscillations.  Cousteix  et  al .  (1977)  used  hot-wire  anemometers  to 
measure  time  averaged  and  phase  averaged  velocity  profiles  as  well  as 
turbulent  intensity  and  shear  stress  across  the  boundary  layer.  They 
found  that  the  general  behavior  of  the  boundary  layer  and  that  the 
structure  of  the  turbulence  were  not  fundamentally  affected  by  the 
unsteadiness  of  the  flow.  A  quasi-steady  approximation,  combined  with 
a  mixing  length  model,  was  used  to  describe  unsteady  flow  at  small 
Strouhal  number.  Later,  Cousteix  et_  al.  (1981,1982)  extended  the 

* 

range  of  their  experiments  so  that  they  covered  Strouhal  numbers  from 
1.5  to  18.  They  confirmed  that  the  mean  flow  field  is  not  affected 
by  the  unsteady  effects,  even  up  to  a  Strouhal  number  18,  An  interesting 
finding  is  that  the  main  parameters  of  boundary  layer  like  displacement 
and  momentum  thickness  6^,  0  ...  oscilate  as  Strouhal  number  increases. 
Direct  measurements  of  the  wall  shear  stress  was  done  by  using  a 
modified  hot  wire  gauge.  They  found  that  the  law  of  the  wall  is  not 
valid  for  evaluating  wall  shear  stress  if  the  Strouhal  number  is  larger 
than  5. 

Parikh  et  al.  (1981,1982)  and  Jayaraman  et  al.  (1982)  designed  a 
test  section  in  which  the  boundary  layer  was  subjected  to  a  sinusoidally 
varying,  adverse,  free-stream  gradient.  In  the  upstream  of  the  test 
section,  a  trubulent  boundary  layer  grew  under  steady,  zero-pressure 
gradient  conditions.  This  provided  a  well-defined  inlet  condition  to 
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the  test  section,  which  is  especially  important  for  computation  of  the 
flow.  The  imposed  oscillation  frequencies  used  in  their  experiments 
covered  a  range  from  quasi-steady  (f  =  0)  to  values  approaching  the 
bursting  frequency.  The  relative  amplitudes  of  oscillations  were  5% 
and  25%  of  mean  free  stream  velocity  at  the  end  of  the  test  section. 

The  laser  Doppler  velocimeter  was  used  to  measure  the  phase  averaged 
velocity  profiles  and  turbulent  intensities.  Again,  they  found  the 
time  averaged  flow  field  is  essentially  unaffected  by  the  imposed 
unsteadiness  and  considered  this  a  consequence  of  two  observations: 

(a)  The  time-averaged  Reynolds  stress  distribution  across  the  boundary 
layer  is  unaffected;  and  (b)  the  Reynolds  stresses  arising  from  the 
organized  velocity  fluctuations  under  imposed  oscillation  are  negligible 
compared  to  the  Reynolds  stresses  due  to  the  random  fluctuations.  How¬ 
ever,  the  periodic  flow  was  observed  to  be  very  strongly  dependent  on 
the  frequency  of  the  excitation.  At  high  frequencies  of  oscillation, 
the  boundary  layer  thickness  and  turbulent  Reynolds  stress  were  found 
to  become  frozen  over  the  oscillation  cycle  while  they  were  observed 
to  have  significant  variation  over  the  cycle  at  low  frequencies.  The 
steady  form  of  the  law  of  the  wall  was  found  not  to  be  valid  for  the 
phase  averaged  velocity  profile  at  high  frequencies.  The  authors 
argued  that  the  Strouhal  number  based  on  the  streamwise  distance  in  the 
test  section  and  the  local  mean  free-stream  velocity  is  a  very  important 
parameter  for  the  periodic  velocity. 

The  unsteady  separating  turbulent  boundary  layer  over  a  flat  plate 
has  been  studied  extensively  by  Simpson  et_  al .  (1981, 1983a, 1983b) .  The 
hot-wire  anemometer,  the  laser  Doppler  velocimeter  and  a  thermal  flow 
direction  probe  were  used  to  measure  the  effects  of  the  sinusoidal 
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unsteadiness  of  the  free-stream  velocity  on  the  separating  turbulent 
boundary  layer  at  low  oscillation  frequency. 

The  measurements  showed  that,  upstream  of  any  flow  reversal  or  back- 
flow,  the  flow  behaves  in  a  quasi-steady  manner.  The  semi-logarithmic 
law-of-the-wall  velocity  profile  was  found  to  be  valid  at  each  phase  of 
the  cycle.  The  phase  and  the  amplitude  of  phase  averaged  velocity  was 
found  close  to  the  values  in  the  free— stream.  However,  downstream  of 
the  detachment,  large  variations  of  the  amplitude  and  phase  of  the  phase 
averaged  velocity  and  the  phase  averaged  turbulent  properties  were 
observed  across  the  boundary  layer.  They  considered  that  unsteady  effects 
produce  a  hysteresis  in  relationships  among  flow  parameters.  The  effect 
of  unsteadiness  was  described  by  the  authors  in  the  following  way:  As 
the  free-stream  velocity,  during  a  cycle,  begins  to  increase,  the  detached 
shear  layer  decreases  in  thickness,  and  the  fraction  of  time  that  the  flow 
moves  downstream  r  increases  as  backflow  fluid  is  washed  downstream.  As 
the  free-stream  velocity  nears  the  maximum  value  in  a  cycle,  the  increas¬ 
ingly  adverse  pressure  gradient  causes  progressingly  greater  near  wall 
backflow  at  downstream  locations  while  r  remains  high  at  the  upstream 
part  of  the  detached  flow.  After  the  free-stream  velocity  begins  to 
decelerate,  the  detached  shear  layer  grows  in  thickness  and  the  location 
where  flow  reversal  begins  moves  upstream.  Higher  oscillation  frequencies 
were  found  to  enhance  this  hysteresis  phenomenon. 

2.3  Periodic  Trubulent  Pipe  (or  Channel)  Flow 

Lu  (1973)  studied  the  effects  of  pulsations  on  a  fully  developed 
turbulent  pipe  flow  of  water  in  a  range  of  Reynolds  numbers  from  16,000 
to  81,600,  with  amplitudes  from  0  to  0.17  of  mean  flow  rate  and  frequencies 
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from  0  to  1.5  Hz.  The  hot-film  anemometer  was  used  to  measure  the 
velocity  profiles  and  turbulent  intensities.  It  was  found  that,  when 
the  pulsating  amplitude  is  small  compared  to  the  mean  velocity  the 
Reynolds  number  and  the  dimensionless  pulsating  frequency  are  two 
important  parameters  to  characterize  the  flow,  whereas,  when  the  amplitude 
is  not  small,  the  pulsating  amplitude  becomes  an  additional  parameter. 

The  time  averaged  flow  field  was  found  not  be  affected  by  the  unsteadi¬ 
ness.  They  suggested  that  a  quasi-steady  model  can  be  used  to  predict 
periodic  flow  at  low  oscillation  frequency. 

Achary  and  Reynolds  (1975)  made  measurements  of  periodic  turbulent 
air  flow  in  a  channel  by  using  a  hot-wire  anemometer.  The  Reynolds  number, 
based  on  channel  half  width,  was  13,800;  the  amplitudes  were  2.4%  and  3.6% 
of  the  channel  centerline  mean  velocity;  the  oscillation  frequencies 
were  24  Hz  and  40  Hz  respectively.  The  data  show  that  the  periodic  velocity 
is  constant  over  most  of  the  flow,  with  all  the  intersting  effects  occur¬ 
ring  in  a  region  very  close  to  the  wall.  The  turbulent  Reynolds  stress 
perturbations  were  found  to  be  out  of  phase  with  the  strain-rates, 
indicating  a  viscoelastic  type  of  repsonse.  An  interesting  finding  is 
the  strange  difference  between  the  25  Hz  data  and  40  Hz  data.  They 
suggested  that  there  might  be  an  interaction  between  the  turbulence  and 
the  imposed  oscillation  since  the  data  40  Hz  was  taken  at  a  frequency  close 
to  the  "bursting  frequency".  Several  closure  models  were  developed  by  the 
authors  for  the  purpose  of  predicting  the  experimental  results.  The 
models  were  not  satisfactory.  However,  they  revealed  several  key  features 
which  need  to  be  incorporated  in  a  successful  model. 

Misushina  et_  al.  (1973)  made  comprehensive  experimental  studies  of 
the  velocity  profiles,  turbulent  intensities  and  auto-correlations  in  a 


14 


2  cm  turbulent  pipe  flow  with  imposed  oscillation  frequencies  from 

4 

0.125  Hz  to  1.32  Hz  and  a  time  averaged  Reynolds  number  10  .  Electro¬ 
chemical  probes  were  used  to  measure  the  instantaneous  velocity  and 
the  turbulent  intensity.  The  bursting  period  was  measured  by  determining 
the  delay  time  of  the  location  of  the  maximum  in  the  auto-correlation 
curve  of  velocity  fluctuations.  The  authors  suggested  that  turbulent 
pulsating  flows  can  be  classified  into  two  groups  by  means  of  a  critical 
period  of  pulsation,  which  is  equal  to  the  maximum  time  between  bursts. 
When  the  periods  of  pulsation  are  longer  than  the  critical  period,  the 
velocity  profiles  are  similar  to  those  for  steady  flow  and  the  intensity 
of  turbulent  fluctuations  don’t  vary  in  a  cycle.  For  shorter  periods, 
the  velocity  profiles  are  different  from  those  for  steady  flow  and  the 
intensity  of  turbulent  fluctuations  oscillates  in  a  cycle  with  its  phase 
about  180°  different  from  that  of  the  phase  averaged  velocity.  Later 
Mizushina  et_  al_.  (1975)  investigated  the  effects  of  pulsations  on  the 
bursting  phenomenon.  They  claimed  that  if  the  pulsating  period  is  in 
the  range  of  the  burst  period,  turbulence  is  generated  near  the  wall  by 
the  flow  pulsation.  Then,  the  generated  turbulence  propagates  radially 
to  the  centre-line  of  the  tube.  The  mean  propagation  time  of  the 
generated  turbulence  from  the  origin  to  the  centre-line  of  the  tube 
agrees  well  with  the  mean  burst  period  of  steady  turbulent  flow  and 
is  independent  of  the  pulsation  period.  In  the  experiment,  they  found 
that  when  pulsation  period  is  less  than  the  mean  burst  period,  the 
turbulent  intensity  oscillates  near  the  wall  and  is  frozen  to  its  steady 
state  value  in  the  central  region  of  the  pipe.  This  was  interpreted 
by  the  authors  to  mean  that  the  oscillation  induced  intensity  near  the 
wall  doesn’t  have  time  to  penetrate  to  the  center  in  a  period. 


Ohmi  £t  al.  (1976)  reported  experimental  results  on  the  velocity 
distribution  in  pulsating  turbulent  pipe  flow.  They  found  that  the 
time-averaged  distribution  in  pulsating  flow  agrees  well  with  that  of 
steady  flows,  and  that  the  oscillating  velocity  profiles  can  be  described 
by  approximate  solutions  derived  from  a  four-region  model  for  a  steady 
turbulent  flow.  They  further  (1978)  solved  the  fundamental  equation 
for  pulsating  turbulent  flow  numerically  by  using  this  four-region  eddy 
viscosity  model.  The  calculation  indicates  that  a  quasi— steady  treatment, 
i.e.,  a  time  dependent  friction  velocity,  should  be  applied  at  low  oscil¬ 
lation  frequencies.  The  numerical  results  also  suggest  that  the  effect 
of  large  Reynolds  number  could  be  similar  to  that  of  low  frequency  on 
the  oscillating  flow  field. 

Kirme  (1979)  used  the  laser  Doppler  velocimeter  to  measure  velocities 
in  pulsating  turbulent  pipe  flow.  He  pointed  out  that  there  might  be  a 
small  region  of  backflow  near  the  wall  for  some  combinations  of  flow 
parameters.  For  a  larger  Reynolds  number,  higher  frequencies  and  amplitud 
of  oscillation  are  needed  to  form  these  zones.  Kirme  also  reported  that 
an  eddy  viscosity  model  based  on  turbulent  kinetic  energy  can  not  predict 
the  pulsating  flow  correctly. 

Kita  at  al.  (1980)  proposed  a  fluctuating  eddy  viscosity  model  for 

Reynolds  stress,  based  on  a  five  region  model  in  steady  turbulent  flow, 

and  compared  it  with  experimental  data.  They  suggested  that  a  modified 

/ci)  * 

Strouhal  number  fi/  /  Re*  ,  where  fl  =  /  v  r„  and  R*  =  ,  characterizes 

U  e  v 

the  volumetric  flow  rate,  while  is  an  important  parameter  near  the 

wall.  Their  experimental  data  show  no  effect  of  oscillation  on  the 
time-mean  velocity  profile. 
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Binder  and  Kuney  (1982)  measured  periodic  velocity  oscillations  near 

the  wall  in  a  turbulent  pulsating  channel  water  flow.  They  reported 

that  with  a  proper  optical  arrangement  it  was  possible  to  use  LDV  to 

make  measurements  as  close  as  0.25  mm  from  the  wall  corresponding  to 

y  =  2  for  u  =17.5  cm/s.  The  results  showed  that  the  mean  flow  and 
c 

the  mean  turbulent  intensity  were  not  affected  by  the  forced  oscillation. 

The  amplitude  and  the  phase  shift  profiles  of  the  periodic  component 

approximately  follow  the  laminar  Stokes  solution  at  high  frequency.  At 

low  frequencies,  the  phase  shift  near  the  wall  decreases  to  slightly 

negative  values  while  the  amplitude  remains  close  to  the  Stokes  solution. 

The  phase  averaged  longitudinal  trubulent  intensity  was  not  found  to  be 

simply  proportional  to  the  velocity  oscillations.  The  authors  argued  that 

in  the  region  near  the  wall  the  Stokes  thickness  nondimensionalized 

s 

with  the  mean  viscous  sublayer  thickness  is  an  important  parameter 
needed  to  characterize  turbulent  pulsating  flow. 

Ramaprian  and  Tu  (1982 ,1983a, 1983b)  reported  detailed  experimental 
data  on  the  instantaneous  velocity  and  wall  shear  stress  at  a  mean 
Reynolds  number  50,000  in  a  fully  developed  turbulent  pipe  flow  in  which 
the  volumetric  flow  rate  was  varied  sinusoidally  with  frequencies  0.5  Hz 
and  3.6  Hz,  and  amplitudes  64%  and  15%  of  the  mean  flow  rate  respectively. 
In  contrast  to  most  other  researchers,  they  concluded  from  their  data 
that  both  the  time  mean  velocity  profiles  and  the  time  mean  turbulent 
intensities  are  affected  by  the  imposed  unsteadiness.  They  considered 
that  these  effects  are  caused  by  the  combination  of  high  oscillation 
frequency  (approaching  the  bursting  frequency)  and  large  amplitude.  The 
phase  averaged  turbulent  intensity  data  showed  that  at  high  frequency 
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the  turbulent  intensity  remains  frozen  throughout  the  cycle  in  the  outer 
region  and  that  at  lower  frequency  it  experiences  significant  variation. 
The  authors  also  reported  the  results  of  direct  measurement  of  wall 
shear  stress  in  turbulent  pulsating  flow  with  flush  mounted  hot  film 
gauges.  The  phase  shifts  of  the  wall  shear  stress  were  found  to  be  of 
the  order  of  10  degrees  in  their  experimental  range.  However,  since 
they  didnf  t  consider  the  frequency  response  of  the  probe,  serious  errors 
could  be  present  in  the  results  obtained  at  high  frequency.  They 
suggested  a  modified  Strouhal  number  ,  which  uses  friction  velocity 
instead  of  bulk  mean  velocity,  to  characterize  unsteady  turbulent  flow. 
With  this  parameter,  they  classified  turbulent  pulsating  flow  into 
five  regimes,  based  on  the  histogram  of  the  intervals  between  the 
bursts  measured  by  Mizushina  et_  al. 

2 . 4  Turbulent  Flow  Over  a  Solid  Wavy  Surface 

A  relevant  field  to  the  pulsating  flow  is  turbulent  flow  over  a 
wavy  surface.  The  spatial  periodic  variation  of  the  solid  boundary 
causes  changes  of  pressure  in  the  fluid  field  and  of  shear  stress  along 
the  wall.  For  small  enough  wave  amplitudes  a  linear  response  is  obtained 
whereby  the  functions  describing  the  pressure  and  shear  stress  variation 
have  a  single  harmonic.  Abrams  et  al.  (1981)  presented  necessary 
conditions  for  a  linear  response. 

A  review  of  the  measurements  of  wall  shear  stress  and  pressure 
profiles  over  a  solid  wavy  surface  is  given  by  Abrams  (1984). 

The  analysis  carried  out  by  Thorsness  et  al.  (1978)  indicates  that 
the  phase  angle  characterizing  the  shear  stress  variation  at  the  wall  is 
a  unique  function  of  a  wave  number,  a  =  2ttv/Au  ,  made  dimensionless 


18 


with  the  kinematic  viscosity,  v,  and  the  friction  velocity,  u  .  The 
amplitude  of  the  shear  stress  variation,  made  dimensionless  by  using 

"k  k 

v  and  u  ,  is  found  to  vary  linearly  with  a^u  /v,  where  a^  is  the  amplitude 
of  the  wave.  The  ratio  of  these  two  dimensionless  amplitudes  is  also  an 
unique  function  of  a+.  To  predict  the  phase  shift  and  amplitude  of  the 
wall  shear  stress,  the  principal  problem  is  the  specification  of  the 
wave  induced  Reynolds  stresses.  Thorsness  (1975)  tried  several  models 
and  found  that  the  mixing  length  model  of  Loyd  et_  al.  (1970)  did  the 
best  job  in  accounting  for  the  influence  of  wave  induced  variations  of 
the  turbulence  properties.  This  model  predicts  a  sharp  variation  of 
phase  shift  of  the  wall  stress  in  a  narrow  range  of  a*\  However,  his 
experimental  data  didnTt  cover  a  large  enough  range  of  conditions  to 
show  this  variation. 

Abrams  (1984)  extended  the  experimental  range  by  using  the  same 

k 

wavelength  X  =  2  inches,  and  by  increasing  the  maximum  value  of  u 

by  a  factor  of  four.  The  experimental  data  indeed  show  a  sharp  change 

-4  +  -3 

in  the  phase  shift  of  the  wall  shear  stress  for  6  x  10  <  a  <  10 

The  interpretation  of  this  observation  is  the  relaxation  effect  of 
turbulent  flow  in  the  viscous  wall  region,  associated  with  the  wave- 
induced  variation  of  the  pressure  gradient.  The  physical  process 
underlying  this  phenomenon  is  not  yet  understood.  Although,  the  author 
tried  a  k-e  model  to  introduce  the  relaxation  in  a  more  natural  way, 
the  calculated  results  couldn’t  describe  the  data  any  better.  His 
work  also  indicates  that  the  wave  induced  variation  of  pressure  gradient 
has  a  more  important  effect  on  the  turbulence  in  the  viscous  sublayer 
than  the  wave  induced  variations  of  streamline  curvature. 
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2.5  Summary  of  the  State  of  Art 

Due  to  large  efforts  made  in  the  last  decade,  many  experimental 
and  some  computational  results  have  been  published.  Nevertheless, 
turbulent  pulsating  flow  is  still  not  understood. 

Most  of  the  data  have  shown  that  the  time  mean  flow  field  is 
unaffected  by  the  imposed  oscillation.  Only  the  data  reported  by 
Ramaprian  and  Tu  indicate  a  slight  effect  of  oscillation  on  the  time 
mean  velocity  profile  and  on  the  time  mean  turbulent  intensity.  One 
can  at  least  get  an  impression  that  the  mean  flow  field  is  not 
sensitive  to  the  imposed  unsteadiness. 

Nevertheless,  the  phase  averaged  velocity  field  is  strongly 

dependent  on  the  characteristics  of  the  external  excitation,  especially 

of  the  oscillation  frequency.  The  response  of  turbulent  flow  to  an 

imposed  oscillation  is  qualitatively  similar  to  that  of  laminar  flow. 

At  very  low  frequency,  the  velocity  profile  for  both  flows  at  each 

instant  is  the  same  as  that  in  a  steady  flow  with  the  same  flow  rate. 

At  high  frequency,  the  periodic  velocity  component  for  both  flows  is 

flat  in  the  center  of  the  pipe  and  has  a  rapid  change  near  the  wall. 

However,  the  response  of  these  two  flows  are  quantitatively  very 

different.  For  example,  the  frequency  parameter  SI  =  r^,  which 

v  0 

characterizes  the  laminar  periodic  flow,  is  no  longer  a  suitable 
parameter  for  the  turbulent  case.  Although  several  frequency  para¬ 
meters  have  been  suggested  for  characterizing  turbulent  pulsating 
flow,  no  definite  conclusion  can  be  deduced.  Moreover,  the  amplitude 
of  the  imposed  oscillation  and  the  mean  Reynolds  number  are  other 
important  factors  affecting  the  turbulent  pulsating  flow.  Their  effects, 
especially  when  they  combine  with  the  effect  of  frequency,  are  not  yet 


understood. 
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There  are  some  data  on  the  phase  averaged  intensity  of  turbulent 
fluctuations  reported  in  the  literature.  It  seems  that  the  phase 
averaged  intensity  oscillates  around  the  mean  value  at  low  frequency, 
while  at  high  frequency  it  varies  only  near  the  wall  and  is  frozen  to 
the  mean  value  in  the  region  away  from  the  wall.  The  phase  of  the 
intensity  variation  also  shifts  from  that  of  the  velocity,  indicating 
that  there  is  no  simple  relation  between  the  turbulent  properties  and 
the  rate  of  strain. 

The  effects  of  the  imposed  oscillation  on  the  turbulent  structure 
are  not  clear.  Most  of  the  existing  data  for  the  Reynolds  stress  are 
calculated  from  the  measurement  of  the  velocity  profile  by  using  the 
momentum  balance  equation.  The  large  possible  errors  in  data  obtained 
in  such  manner  make  it  difficult  to  come  to  any  definite  conclusions. 

A  strong  interation  between  the  turbulent  and  the  imposed  oscillations 
at  some  resonant  frequency,  anticipated  by  several  researchers,  is 
still  not  substantiated. 

The  response  of  the  flow  in  near  wall  region  is  very  important  for 
understanding  turbulent  pulsating  flow.  Because  the  inertial  effect 
is  much  less  in  this  region  than  in  the  outer  flow,  the  data  near  the 
wall  offer  a  good  test  for  turbulent  models  used  to  predict  the  flow. 
However,  few  reliable  measurements  have  been  made  in  this  region.  At 
high  frequencies,  the  layer,  in  which  the  periodic  velocity  has  spatial 
variation,  is  so  thin  that  measurement  of  its  properties  is  not  feasible. 
Therefore,  the  wall  shear  stress  measurements  in  the  present  study  are 
needed  to  obtain  a  better  understanding  of  a  turbulent  pulsating  flow. 
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CHAPTER  3 .  THEORY 

In  this  chapter,  a  theoretical  framework  for  the  description  of 
turbulent  pipe  flow  with  a  superimposed  organized  sinusoidal  excitation 
is  presented.  The  flow  field  is  decomposed  into  three  components:  the 
time  mean  flow,  an  oscillation  flow  and  the  turbulent  fluctuations. 

The  equation  governing  the  oscillation  flow  is  derived  and  solved 
numerically  by  applying  several  turbulent  models  based  on  the  eddy 
viscosity  concept.  Then,  the  phase  and  the  amplitude  characterizing 
the  variation  of  shear  stress  at  the  wall  are  calculated  at  different 
flow  conditions.  The  results  are  compared  with  experimental  data  in 
Chapter  5. 

3.1  Triple  Decomposition  and  the  Phase  Average 

In  steady  turbulent  flow,  the  problem  is  attacked  by  decomposing 
flow  variables  into  time  mean  values  and  turbulent  fluctuations.  The 
term  "steady”  means  the  time  mean  values  are  independent  of  time. 

When  an  organized  unsteadiness  is  imposed  on  the  system,  this  additional 
component,  which  is  a  deterministic  function  of  time,  has  to  be  included. 
Following  Hussain  and  Reynolds  (1970) ,  the  problem  can  be  formulated  as 

E(x,t)  =  E(x)  +  E(x,t)  +  E ! (x, t) ,  (3.1) 

where  E(x)  is  the  time  mean  value  of  E(x,t)  at  a  given  location  x, 

E(x, t),  the  organized  response  component  due  to  the  imposed  excitation 
and  ET(x, t),  the  turbulent  fluctuation. 

In  order  to  separate  the  three  components  of  a  flow  variable,  two 
different  averaging  procedures  are  employed;  i.e.,  a  time  average  and  an 
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ensemble  average.  The  time  average  is  conventionally  used  in  steady 
turbulent  studies,  and  is  defined  as 

1  fT 

E(x)  =  lim  —  E(x, t)  dt  (3.2) 

T-> 00  1  J0 

In  discrete  form,  the  time  average  is  the  average  of  a  large  number 
of  samples  in  a  record  of  sufficiently  long  duration.  That  is 

1  N 

E(x)  =  lim  -  £  E  (x,t1)  (3.3) 

N  ^  oo  i=l 

In  practice,  the  total  time  duration  in  which  all  the  samples  were  taken 
is  large  compared  to  the  imposed  excitation  period  and  the  longest  turbulent 
fluctuation  period. 

The  ensemble  average  is  defined  as  the  average  over  a  large  number 
of  realizations  of  the  same  phenomenon  and  is  expressed  as 

1  N 

<E(x,t)>  =  lim  -  l  Ei(x,t)  ,  (3.4) 

N->  0°  i=l 

where  i  identifies  a  particular  realization.  This  procedure  can  be 
used  to  extract  the  organized  flow  component. 

When  the  imposed  excitation  is  periodic,  the  ensemble  average  is 
equivalent  to  a  phase  average.  Each  cycle  is  considered  as  a  realization 
in  the  ensemble  average.  So,  a  phase  average  is  defined  as  the  average 
of  samples  from  a  large  number  of  cycles  at  a  fixed  cycle  angle. 

N-l 

<E(x,0)>  =  lim  —  £  E(x,0  +  27rn ) ,  0  <  0  <  2tt, 

N->  °o  n=Q 


(3.5) 
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where  0  is  the  cycle  phase  angle.  In  practice,  N  is  a  large  finite 
number. 

By  definition,  phase  averaging  eliminates  the  term  involving  random 
turbulent  fluctuations,  namely, 

<E(x,0)  >  =  E(x)  +  E(x,0)  ,  (3.6) 

and  the  time  average  of  <  E(x,0)  >  eliminates  the  oscillation  component, 

<  E(x,0)  >  =  E(x)  .  (3.7) 

From  the  above  expressions,  the  oscillation  component  and  the 
turbulent  fluctuations  can  be  extracted  from  the  flow  variable  as 


E(x,t)  =  <  E(x,t)  >  -  E(x)  ,  (3.8) 

and 

E'(x,t)  =  E(x,t)  -  <  E(x, t)  >  .  (3.9) 


Figure  3.1  displays  the  relations  among  these  components. 


3.2  Governing  Equation  for  the  Oscillation  Component 


The  governing  equations  for  an  incompressible  flow  of  a  Newtonian 
fluid  in  a  circular  pipe  are  given  by  the  Nevier-Stokes  equation  in 
cylindrical  coordinates, 
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Figure  3.1  Time  and  phase  averages 
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along  with  the  continuity  equation 


19  ,  n,19w,9w 

rJi  (rv)  +r  36  +  3x  =  0 


(3.12) 


(3.13) 


where  u,  v,  and  w  are  the  velocity  components  in  the  axial  direction  x, 
the  radial  direction  r  and  the  angular  direction  0  respectively. 

For  turbulent  pulsating  flow,  the  triple  decomposition  is  applied 
to  these  flow  variables  as, 

u  =  u  4*  u  +  u’  , 

v  =  v  +  v  +  v*  , 

w  =  w  4*  w  4*  wf  , 

p  =  p  +  p  +  p'  .  (3.14) 

If  the  flow  is  fully  developed  and  symmetric  in  the  angular  direction, 
time  averages  of  equations  (3.10)  to  (3.14)  are  dramatically  simplified. 

The  components  of  velocities  such  as  v,  v,  w,  w  vanish.  The  equation 
describing  the  time  mean  flow  becomes 


0  = 


_  _ 

dx 


1  _d 
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+  -f 

r  dr 


(3.15) 


To  solve  this  equation,  the  Reynolds  stress,  x  =  -pu'v',  has  to 
be  specified  in  terms  of  u  and  r.  This  is  the  so  called  closure  problem. 
In  turbulent  pulsating  flow,  one  principal  question  is  whether  the  time 
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mean  turbulent  property  u'v'  is  affected  by  the  imposed  unsteadiness. 

If  it  is  not  affected,  the  closure  method  used  for  steady  flow  can  be 
applied  and  the  time  mean  velocity  profile  will  be  the  same  as  that  for 
steady  flow.  Otherwise,  the  time  mean  velocity  profile  in  pulsating 
flows  should  show  some  differences  from  that  in  steady  flow,  and  a  new 
closure  method  needs  to  be  developed. 

In  a  similar  way,  the  equation  describing  phase  averaged  flow  is 
derived  by  phase  averaging  equations  (3.10)  to  (3.13).  It  is  given  as 


3<  u  > 

at 


i±2^  _  I  _L  (r  <  T(t)>)  +  £  JL 

3x  r  3r  '  +  r  3r 


r3<  u  > 


3r 


(3.16) 


where  <r  >  = -p  <  u'v'  >  is  the  phase  averaged  Reynolds  stress. 

The  equation  which  governs  the  oscillation  component  is  obtained 
by  subtracting  the  mean  flow  equation  (3.15)  from  the  phase  averaged 
equation  of  motion  (3.16).  This  gives 
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Xu  order  to  solve  equation  (3.17)  the  oscillation  induced  Reynolds  stress 
t^,  which  is  defined  as 


i<‘>  .  <T<‘>>  -7(t> 


(3.18) 


has  to  be  specified.  This  is  done  by  following  Boussinesq’s  (1877) 
concept  of  turbulent  eddy  viscosity,  which  relates  Reynolds  stress  to 
the  velocity  gradient  by  an  analogy  with  Newton’s  law  of  viscosity.  The 
phase  averaged  Reynolds  stress  then  can  be  written  as 
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If  it  is  assumed  that  the  imposed  small  oscillation  induces  a  variation 
of  turbulent  eddy  viscosity  v  around  its  time  mean  value  vfc,  the 
oscillation  induced  Reynolds  stress  is  expressed  as 
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3u 
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(3.20) 


if  it  is  assumed  that  the  second  order  term  v  7—  is  negligible.  The 

t  <3y 

time  mean  Reynold  stress  is  written  as 


7(t)=  P 
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(3.21) 


3.3  Time  Mean  Flow 

Most  previous  researchers  have  found  that  the  time  mean  flow  is  not 
affected  by  imposed  small  amplitude  flow  oscillations.  Experimental 
results  on  wall  shear  stress  obtained  in  the  present  study  also  support 
this  conclusion.  Therefore,  the  time  mean  velocity  profile  is  calculated 
by  using  the  same  eddy  viscosity  v  that  represents  a  steady  flow.  That 
is, 
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(3.22) 


where  u  is  the  time  mean  friction  velocity  and  Tq  is  the  radius  of  the 
pipe.  An  empirical  equation  proposed  by  Reichardt  (1951)  is  used  for 
the  prediction  of  v  in  the  core  region: 
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where  k  is  the  von  Karman  constant,  0.4. 


28 


In  the  region  close  to  the  wall,  the  mixing-length  theory  is  used 
so  that 
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t 


5u 

ay  ’ 


(3.24) 


where  £.,  the  mixing  length  characterizing  the  undisturbed  flow,  is 
described  by  the  Van  Driest  function, 


£  =  ky  [1  -  exp(-D  )] 

m 


(3.25) 


The  term  is  a  damping  function  defined  as 
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and  is  the  time  mean  shear  stress  at  the  wall.  The  Van  Driest  damping 
constant  A  is  26  for  pipe  flow.  By  combining  (3.24)  and  (3.25),  the 
eddy  viscosity  in  the  wall  region  can  be  expressed  as 
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(3.27) 


Figure  3.2  gives  the  radial  variation  of  eddy  viscosity,  calculated 
from  (3.23)  and  (3.27).  It  is  seen  that  the  change  of  the  eddy  viscosity 
from  the  near  wall  region  to  the  outer  flow  is  fairly  smooth.  The 
velocity  profile  calculated  from  (3.22)  at  Re  =  23,000  is  compared  with 
the  data  measured  by  Nikuradse  (1932)  in  Figure  3.3.  The  good  agreement 
between  the  two  indicates  that  the  eddy  viscosity  model  used  in  this 
study  is  able  to  describe  the  time  mean  flow  field.  Furthermore,  the 
Van  Driest  formula  is  a  continuous  function  from  the  wall  through  the  viscous 
sublayer  and  into  the  fully  developed  turbulent  region.  This  feature 


EDDY  VISCOSITY  (v  /v) 


29 


Y/R 

Figure  3.2  Turbulent  eddy  viscosity  distribution  for 
steady  flows  (calculated  from  (3.23)  and 
(3.27)) 
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Figure  3.3  Time  mean  velocity  profile  of  steady 

turbulent  flow  [Re  =  23,000,  data  from 
Nikuradse  (1932) ,  —  calculated  from  (3.22)] 
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is  very  important  in  modeling  an  oscillating  flow  for  which  the  flow 
variation  at  high  frequency  case  is  confined  to  a  thin  layer  near  the 
wall. 


3.4  Quasi-laminar  Model  (Model  A)  for  the  Induced  Flow  Oscillation 


The  simplest  model  for  the  induced  flow  oscillations  involves  the 
assumption  that  induced  oscillations  of  the  Reynolds  stress  are  zero, 

"r  =  0.  This  means  that  the  induced  oscillations  in  a  turbulent 
flow  will  behave  as  the  same  as  for  a  laminar  pulsating  flow,  i.e.. 


3u  _  _  .1  j}jp  v_  _3_  f  3u 

9t  p  3x  r  3r  (r  3r  j 


(3.28) 


This  equation  has  been  solved  by  Th.  Sexl  (1930)  and  S.  Uchida  (1956) 
for  a  pressure  gradient  with  sinusoidal  variation  in  time. 


1  3p 

- t—  =  a  cos  ojt  . 

p  3x 


(3.29) 


The  model  is  assumed  to  approximate  the  flow  at  very  high  frequencies, 
for  the  spatial  variation  of  the  imposed  flow  oscillations  is  confined 
to  a  very  thin  layer  near  the  wall  where  turbulence  has  negligible 
effect.  At  high  frequencies,  ft  =  rQ  »  1,  the  solution  is  given  as 


u(r , t) 


a)  • 

—  <  sin  wt 
“1 


sin 


wt  - 


/  W  . 

/  2V  (r 


0 


(3.30) 


This  equation  shows  that  the  oscillation  component  in  the  central 


portion  of  the  pipe  behaves  as  a  solid  body,  i.e., 
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u(t) 


r-K) 


=  —  sin  cot 
co 


(3.31) 


Only  in  the  Stokes  layer,  of  thickness 


(3.32) 


does  the  velocity  u  (r,t)  change  rapidly  with  distance  y  from  the  wall. 
The  velocity  gradient  at  the  wall  is  obtained  from  (3.30)  as 


3u 

3y 


y=0 


a 

—  cos 

COV 


(3.33) 


Certain  phase  relations  are  observed  from  these  equations:  The  central 
velocity  lags  the  pressure  gradient  by  90°.  The  velocity  gradient  at 
the  wall  lags  by  only  45°,  i.e.,  the  velocity  gradient  at  the  wall  is 
45°  ahead  of  the  central  velocity. 

Figure  3.4  shows  the  oscillation  velocity  profile  at  eight  phases 
for  laminar  flow  when  the  dimensionless  frequency  parameter  ft  =  150. 

The  oscillation  velocity  component  is  made  dimensionless  with  the  amplitude 
of  velocity  at  the  center  line.  It  is  seen  that  the  velocity  is  flat  at 
Y  >  0.04  for  each  instant.  The  changes  of  the  velocity  happen  in  a  thin 
layer  Y  <  0.04.  The  crossover  of  the  velocity  profiles  at  different 
phases  in  the  near  wall  region  indicates  the  phase  shift  of  the  velocity. 
The  frequency  effects  on  the  oscillation  velocity  can  be  easily  seen 
in  Figure  3.5  and  Figure  3.6,  where  the  amplitude  and  the  phase  of  the 
oscillation  component  are  plotted  for  four  different  ft.  As  ft  increases, 
the  variation  of  the  oscillation  component  is  confined  to  a  thinner 
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Y*10 


Figure  3.4  Oscillation  velocity  profile  at  different 
phases  for  laminar  flow  (53  =  150) 


AMPLITUDE 


Figure  3.5  Amplitude  profile  of  laminar  oscillation 
flow  at  different  ft 
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Figure  3.6  Phase  profile  of  laminar  oscillation 
flow  at  different  ft 
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region  close  to  the  wall.  An  interesting  result  is  that  the  phase  of 
the  velocity  near  the  wall  is  ahead  of  that  in  the  center. 

A  physical  interpretation  of  the  behavior  of  the  oscillation  velocity 
in  a  laminar  flow  is  given  by  Lighthill  (1954)  and  Gerrard  (1971).  As 
the  piston  oscillates,  the  variation  of  pressure  gradient  is  imposed 
on  the  flow.  Due  to  inertia,  the  change  of  velocity  lags  behind  that 
of  pressure  gradient.  However,  the  inner  part  of  the  flow  responds  to 
the  pressure  gradient  more  rapidly  than  the  outer  flow  since  the  inertia 
of  the  fluid  is  less  here  than  at  the  pipe  center.  The  vorticity, 
which  is  simply  the  radial  gradient  of  velocity  in  a  fully  developed 
pipe  flow,  created  at  the  wall  diffuses  into  the  flow.  The  relevant 
parameter  ft  is  a  measure  of  the  ratio  of  the  tube  radius  to  the  distance 
vorticity  diffuses  in  one  period  of  oscillation.  At  high  frequency  the 
vorticity  never  has  time  to  diffuse  to  the  center  region  before  being 
annulled  by  oppositely  signed  vorticity.  This  results  in  the  flat 
distribution  in  the  center  part  of  the  pipe.  On  the  other  hand,  the 
sectional  mean  velocity  must  equal  to  the  piston  speed  at  each  instant 
by  the  continuity  condition  of  incompressible  fluid.  Retardation  due 
to  viscous  effects  near  the  wall  is  balanced  by  an  acceleration  over 
the  rest  of  the  cross  section.  This  causes  an  overshoot  of  velocity 
at  the  outer  part  of  the  Stokes  layer. 

3.5  Hussain  and  Reynolds  Model  (Model  B)  for  the  Induced  Flow  Oscillations 

If  the  effect  of  oscillation  on  turbulence  is  considered,  i.e., 

f  0,  a  simple  way  to  relate  the  oscillation  induced  Reynolds  stress 
to  the  induced  velocity  gradient  is  to  assume  that  the  eddy  viscosity  is 
specified  only  by  mean  flow  field  (Hussain  and  Reynolds,  1970).  That  is 
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~(t)  —  9u 

T  =  pVt  87 


(3.34) 


The  substitution  of  (3.34)  into  (3.17)  gives 


9u  _  1_  v  _9_ 

9t  p  9x  r  9r 


r  r 


l1 


3u 

3r 


(3.35) 


where  v  is  calculated  by  (3.23)  and  (3.27)  for  core  and  wall  region, 
respectively. 

3.6  Quasisteady  Approximation  for  the  Mixing  Length  (Model  C) 

The  shortcoming  of  Model  B  is  that  it  doesn’t  take  account  of  the 
changes  of  eddy  viscosity  during  the  oscillation.  To  consider  this, 
the  mixing  length  is  assumed  to  be  given  by 

i 

£  -T  +  i  ,  (3.36) 

where  l  is  the  mixing  length  for  the  undisturbed  flow  and  i,  the 
oscillation  induced  part.  Then,  the  oscillation  induced  eddy  viscosity 
can  be  expressed  as 


v  =  i2  -jp  +  211  |H. 
t  9y  9y 


(3.37) 


where  l  needs  to  be  specified. 

The  Van  Driest  formula  for  the  mixing  length  in  (3.25)  is  for  a 

steady  flow  with  zero  pressure  gradient.  The  wall  shear  stress  ~  is 

w 

used  in  the  damping  factor  since  the  shear  stress  is  nearly  constant  in 
the  region  close  to  the  wall.  In  flow  over  a  wavy  surface,  the  wave- 
induced  pressure  gradient  can  cause  large  variations  of  the  shear  stress 
through  the  viscous  wall  region  which  alter  the  production  of  turbulence. 
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Abrams  and  Hanratty  (1984)  took  this  into  account  by  redefining  the 
damping  factor  as 


1  2, 

D  =  y_i _ (z I 

M  1/2  . 

p  vA 


(3.38) 


where  the  local  shear  stress,  x(y),  is  used  instead  of  the  value  at  the 

surface,  x  .  This  formulation  will  predict  a  thickening  of  the  viscous 

wall  region  in  favorable  pressure  gradients  since  x(y)  is,  then,  less 

than  the  value  at  the  surface  x  . 

w 

In  a  pulsating  flow,  the  time  variation  of  the  pressure  gradient 
also  causes  a  large  variation  of  shear  stress  in  the  viscous  wall  region. 
Therefore  local  phase  averaged  shear  stresses  <x(y,t)  >  are  used  in  the 
damping  factor.  The  oscillation  induced  mixing  length  l  is  then  given 
as 


i  = 


y+  t (y , t) 

A  7(y) 


exp 


(-  y+/A) , 


(3.39) 


where 


T(y,t) 


~  (t)  /  .  9u 

"  T  (y’c)  +  y 


(3.40) 


and 


t (y)  =  xw  (1  -  y/rQ) 


(3.41) 


From  (3.20),  (3.37),  and  (3.39)  to  (3.41),  the  oscillation  induced  eddy 
viscosity  is  expressed  as 


v 


t 


vt  +  B  (v  +  vt) 


(1  -  B) 


3u 

3y 


3u 

3y 


(3.42) 


where 


2 


(3.43) 


_  + 
B  =  ky l 


exp  (-  y+/A)/[(l  -  y/r0)  u  *2  ] 


is  a  function  of  y,  has  a  maximum  near  the  wall  and  approaches  to  zero 
with  increases  of  distance  from  the  wall,  as  shown  in  Figure  3.7.  The 
oscillation  induced  Reynolds  stress,  as  described  in  (3.20),  can  be 
written  as 


l(t)  .  ,  x  3u 

T  =  p  h(y)  37 


(3.44) 


with 


^(y)  = 


2vt  +  Bv 
1  -  B 


(3.45) 


This  model  is  a  quasi-steady  approximation.  The  mixing  length, 
as  well  as  the  Reynolds  stress ,  are  directly  associated  with  the 
instantaneous  vleocity  gradient  in  the  same  way  as  for  a  steady  flow. 

3-7  Relaxation  Consideration  of  Pressure  Gradient  Effect  (Model  D) 

Loyd,  Moffat  and  Kays  (1970)  and  Reynolds  (1974),  (1976)  found  that 

(3.26)  with  A  =  25  does  not  correctly  predict  the  effect  of  pressure 

gradient  in  equilibrium  boundary  layers.  They  therefore  allowed  A  to 

be  a  function  of  the  dimensionless  pressure  gradient  p+  =  ^  v/ou*3 

dx 

For  small  p  , 

A  =  A  (1  +  k±  p+  +  k2p+2  +  ...)  .  (3.46) 

with  k^  =  -30  and  k2  =  1.54  x  103. 

For  situations  in  which  the  pressure  gradient  is  varying  rapidly, 

Loyd  et  al. ,  (1970)  suggest  that  the  flow  close  to  the  wall  sees  an 
effective  pressure  gradient  given  by  the  equation 


B  (Y) 
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Figure  3.7  Profile  of  B(Y) 


(Re  =  30,000) 
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r  > 
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[dxj 
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d  (xu*/v)  k 

-Li 


(3.47) 


where  is  a  relaxation  constant  approximately  equal  to  3000. 

Thorsness  et  al.  (1978)  and  Abrams  (1984)  adopted  this  approach  in 
accounting  for  the  non- equilibrium  effects  in  the  viscous  wall  region 
for  wavy  surface  flows.  With  =  -35,  kL  =  1800,  Abrams  obtained 
good  agreement  between  the  model  and  experimental  measurements  of  the 
variation  of  wall  shear  stress  along  a  wave  surface. 

The  pressure  gradient  in  a  pulsating  flow  is  varying  with  time 
rather  than  with  space.  A  straightforward  application  of  (3.47)  to 
this  case  requires  the  definition  of  a  convection  velocity  which 
characterizes  the  streamwise  propagation  of  disturbances  in  viscous 
wall  region  to  relate  spatial  variations  to  time  variations.  A  convection 
velocity  of  C*  =  15  has  been  chosen  to  take  account  of  this  effect.  There- 
fore,  the  effective  pressure  gradient  for  a  pulsating  flow  is  defined  as 
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eff 
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dpi 

dxj  eff 


V15 


(3.48) 


The  oscillation  induced  eddy  viscosity 


is  then  given  as 


v 


t 


Jt  +  B(V  +  V  3u  _  2B  (1  -  y/rQ)  u  *2  £ 

-  B>  f  ^  <1  -  ®>f*  1 


(3.49) 


eff*  A  comparison  with  (3.42)  shows  that  the  second  term 
A 

on  the  right  side  of  (3.49)  considers  the  influence  of  flow  relaxation 
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for  high  frequencies.  In  general,  this  term  is  not  in  phase  with  the 

oscillation  velocity  gradient  .  It  is  a  strong  contributor  to 

in  the  wall  region.  The  oscillation  induced  Reynolds  stress  is  described 


as 


;<*>  = 


*i(y)  |f  -  *2<y)  f 

A  • 


(3.50) 


where 


2B  (1  -  y/r  )  u  2 

My)  = - - - 

1  -  B 


(3.51) 


These  relations  predict  that  the  relaxation  effect  is  confined  in 
the  wall  region  because  B  approaches  zero  at  large  y. 

In  the  core  region,  Z  is  assumed  to  be  zero.,  and  =  2pv^  —  , 

where  v  is  taken  as  a  constant  calculated  from  (3.23)  at  the  pipe 
center.  The  match  point  y^  is  where  the  calculated  ^(Yq)  Is  equal  to 
twice  the  value  of  the  eddy  viscosity  evaluated  at  pipe  center. 


3.8  Numerical  Scheme 


The  time  variation  of  the  wall  shear  stress  caused  by  an  imposed 
variation  of  pressure  gradient  defined  by  (3.29)  is  calculated  by 
solving  equation  (3.17)  numerically.  The  equation  was  made  dimensionless 
in  the  following  manner: 


If  =  cos  T  +  -j 


0(1  -  Y) 
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(3.52) 
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with 


U  = 


Y  = 


T  =  wt. 


(3.53) 


and 

u'  . 

a 

The  boundary  conditions  are 

U  =  0  at  Y  =  0 


and 


^  =  0  at  Y  =  1  (3.54) 

The  equation  was  solved  numerically  by  a  Crank-Nicolson  implicit  scheme 
with  variable  spatial  grid  size.  Since  the  oscillation  velocity  changes 
rapidly  in  the  wall  region  and  the  final  goal  of  the  calculation  is  to 
find  the  velocity  gradient  at  the  wall,  it  is  necessary  to  use  a  very 
small  grid  size  close  to  the  wall.  Grid  sizes  which  are  increasing  in 
a  geometric  series  of  Y  were  chosen: 

AY  (H^_1  -  1) 

Y.  =  -  ,  l*j<J  (3.55) 

J  (H  -  1) 

where  AY^  is  the  first  grid  size  near  the  wall,  and  factor  H  is  the 
ratio  of  two  consequent  grid  sizes.  Ratio  H  is  chosen  as  1.05  for  present 


work.  The  total  number  of  grid  points  in  the  radial  direction  is  then 
given  as 


J  =  log 


+  1 


/  log  H  , 


(3.56) 


where  J  takes  integer  number. 

The  first  grid  size  AY^  was  selected  small  enough  so  that  the 
velocity  gradient  at  the  wall  can  be  calculated  by  the  assumption  of 
linear  variation  of  velocity  in  that  region.  After  several  tries,  AY^, 
was  taken  as  1/30  of  the  Stokes  layer  thickness.  Figure  3*8  gives  the 
schematic  description  of  the  grids. 

The  initial  condition,  time  zero,  was  chosen  as  the  quasi-laminar 

2tt 

result.  Equal  time  steps  of  AT  =  -rrr~  were  used.  At  each  new  time  level 

IZo 

the  velocities  at  all  points  in  Y  direction  were  calculated  from  previous 
time  level  results.  After  about  five  periods  calculation,  the  velocities 
at  each  phase  of  the  cycle  are  found  to  converge  to  constant  values. 

Because  of  the  variable  grid  size  in  space,  the  finite  difference 
expression  for  derivatives  used  with  an  equal  grid  size  cannot  be 
employed.  Expressions  for  finite  difference  with  variable  grid  size 
have  been  derived  by  approximating  a  function  piecewise  with  a  second 
order  polynomial 


U(Y) 


a.  +  b.Y  +  C.Y2, 
3  3  3 


at  Y,  . 
3-1 


(3.57) 


where  U(Y)  is  the  approximated  value;  a^  ,  b_.  and  c_.  ,  the  coefficients 
in  the  polynomial.  Then,  the  finite  difference  forms  of  derivatives 


are  derived  as 
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Figure  3.8  Grids  used  for  finite  difference 


and 
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(3.59) 


where  2  £  j  <  J  -  1.  Details  of  the  derivation  are  given  in  Appendix  II, 
It  is  clear  that  the  above  forms  reduce  to  the  standard  central  finite 


difference  expressions  for  equal  grid  sizes,  i.e.,  H  =  1. 
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CHAPTER  4.  DESCRIPTION  OF  EXPERIMENTS 

4.1  Experimental  Flow  Loop 

A  schematic  drawing  of  the  loop  in  which  the  experiments  were  performed 
is  shown  in  Figure  4.1.  The  loop  was  originally  built  by  Sirkar.  Details 
regarding  its  design  may  be  found  in  a  thesis  by  Sirkar  (1969)  and  in  a 
paper  by  Sirkar  and  Hnaratty  (1970).  In  this  work,  a  reciprocating  piston 
was  connected  to  the  loop. 

The  mean  flow  was  generated  by  a  centrifugal  pump  while  the  imposed 
oscillation  was  introduced  by  a  piston  with  a  diameter  of  6.35  cm  and  an 
adjustable  stroke  length  from  0  to  15.2  cm.  The  mean  flow  rate  was 
measured  by  a  vortex  shedding  flow  transmitter  (Eastech,  Model  2420) 
with  an  accuracy  of  0.5%.  The  frequency  of  oscillation  in  the  experiment 
was  fixed  at  0.325  Hz  or  0.625  Hz  by  adjusting  the  circular  frequency  of 
the  motor  driving  the  piston.  The  oscillating  flow  component  was  intro¬ 
duced  into  the  system  between  diaphragm  valve  3  and  calming  section  4. 

The  uniform  flow  that  emerges  from  the  calming  section  is  tripped  by  a 
1.26  cm  long  ring  consisting  of  a  series  of  0.95  cm  equilateral 
triangles  around  the  circumference  of  the  entrance  of  the  19.4  cm  pipe 
line.  In  the  experiment,  diaphragm  valve  3  and  bypass  valve  17  were 
adjusted  to  impose  a  large  pressure  drop  across  them.  This  minimized 
the  distortions  of  the  sinusoidal  oscillation  in  the  19.4  cm  pipe, 
caused  by  the  time  response  of  the  whole  system. 

Rectangular  electrode  made  of  platinum  were  mounted  flush  with  the 

wall  of  a  test  section  located  67.5  pipe  diameters  downstream  from  the 

pipe  entrance.  The  fluid  was  an  aqueous  solution,  0.1  M  in  potassium 

iodide  and  0.0038  M  in  iodine,  maintained  at  a  temperature  of  25  ±  0.1°C. 

o 

The  kinematic  viscosity  of  the  solution  was  0.00866  cm  /sec. 
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1.  Storage  tank 

2.  Worthington  pump 

3.  Diaphragm  valve 

4.  Calming  section 

5.  19.4  cm  pipe  upflow 

6.  Test  section 

7 .  Anode 

8.  Vortex  flow  meter 

9 .  Temperature 
control  unit 

10.  Recycling  pump 

11.  Filter  pump 

12.  Filter 

13.  10.2  cm  pipe  downflow 

14.  10.2  cm  bypass  pipe 

15.  Pulsating  generator 

16.  Pressure  transducer 

17.  Bypass  valve 


Figure  4.1  Schematic  of  19.4  cm  flow  loop 
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4. 2  Electrochemical  Technique 

Electrochemical  techniques  were  used  to  measure  flow  in  the  immediate 
vicinity  of  the  solid  wall.  These  were  developed  by  Reiss  and  Hanratty 
(1963)  in  this  laboratory.  The  principle  behind  them  is  as  follows: 

When  a  voltage  is  applied  to  an  electrochemical  cell  in  an  aqueous 
solution  of  potassium  iodide  and  iodine,  the  following  reactions  are 
carried  out  on  the  surfaces  of  the  electrodes: 

Cathode  +  2e  +  31 

Anode  31  +  +  2e. 

When  the  electrochemical  process  is  controlled  by  the  rate  of  mass 
transfer  at  the  surface  of  the  cathode  (the  test  electrode) ,  the  current 
flowing  in  the  cell,  I,  is  related  to  the  rate  of  mass  transfer,  N,  by 
the  equation 


N  = 


I 


A  n  F  ’ 
e  e 


(4.1) 


Where  A  is  the  area  of  the  test  electrode,  n  ,  the  number  of  the 
e  e 

electrons  involved  in  the  reaction  and  F,  Faraday Ts  constant. 

Consider  a  rectangular  probe  embedded  in  the  wall  with  its  long 
side  perpendicular  to  the  direction  of  mean  flow.  If  the  dimensions 
of  the  element  are  selected  correctly  (Hanratty  &  Campbell,  1983), 
diffusion  in  both  the  transverse  and  streamwise  directions  may  be 
neglected.  Furthermore,  the  thickness  of  the  concentration  boundary- 
layer  will  be  small  enough  that  the  effect  of  velocities  normal  to  the 
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surface  can  be  neglected  and  the  velocity  in  the  streamwise  direction 
is  given  by 


u  =  Sy 


(4.2) 


The  mass  balance  equation  in  the  concentration  boundary  layer  can  be 
described  as 


3C  9C 

+  Sy  37  = 


32C 

3y2 


(4.3) 


where  a  is  the  diffusivity.  Equation  (4.3)  is  to  be  solved  using  the 
boundary  conditions 


C  (x,  o,  t)  =  0  0  ^  x  *  L 

C  (x,  t)  =  C  (0,  y,  t)  =  Cb  (4.4) 

For  steady  flow  or  if  the  flow  is  unsteady  and  changing  slowly 
with  time,  the  first  term  in  (4.3)  can  be  neglected.  The  solution  of 
this  equation  is  the  pseudosteady  state  approximation  given  by  Reiss 
(1963): 


1 

T(4/3) 


(4.5) 


where 


q  =  y 


+  c  1/3 
y  Sc 


(4.6) 


and  Sc,  the  Schmidt  number.  The  terms  y"*"  and  x+  are  the  coordinate 

k 

distances  made  dimensionless  with  respect  to  friction  velocity  u 
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and  kinematic  viscosity  v.  The  transfer  coefficient,  K  =  N/C,  ,  is 

b 


then  given  as 


(4.7) 


y=0 


Equation  (4.7)  is  normally  used  to  obtain  the  velocity  gradient  at  the 
wall,  or  wall  shear  stress,  from  the  measured  values  of  K.  Details 
regarding  the  measurement  of  the  current  and  the  electrochemical  technique 
are  given  in  a  review  article  by  Hanratty  and  Campbell  (1983). 

When  this  technique  is  applied  to  unsteady  flows,  especially  at 
high  frequencies  of  oscillation,  the  pseudosteady  approximation  might 
not  hold.  Fortuna  and  Hanratty  (1971)  considered  the  effect  of  the 
frequency  response  of  concentration  boundary  layer  on  the  measured 
amplitude  of  the  fluctuating  flow  by  solving  a  linearized  form  of  the 
time  varying  conservation  equation.  In  the  present  work,  particular 
attention  was  given  to  the  phase  lag  errors  that  arise  in  using  (4.7) 
to  calculate  time  variation  of  xw(t) .  Numerical  and  experimental 
studies  of  this  problem,  summarized  in  a  recent  paper  by  Mao  and 
Hanratty  (1984) ,  are  described  in  the  following  section. 

4.3  Use  of  Electrochemical  Wall  Probes  to  Pulsating  Flow 

Let  C  and  S  be  the  sum  of  average  and  time  dependent  terms. 


C  =  C  +  c 


S  =  S  +  s 


(4.8) 


52 


Equation  (4.3)  can  then  be  decomposed  into  equations  for  C  and  for  c. 
The  solution  for  C  is  the  same  as  (4.5).  The  equation  for  c  is  solved 
by  using  a  harmonic  variation  for  s, 


-  iwt 

s  =  s  e  , 

(4.9) 

where  s  and 

is  linear, 

a)  are  real  numbers.  Since  the  differential 

equation  for  c 

a 

^  irnt  s  iajt 

c  =  c  e  =  c  —  e  , 

1  S 

(4.10) 

where  c  and  can  be  complex  numbers. 

The  differential  equation  defining  c^  is 

.  *  ,+2/3  ,  *  3C1  .  *  9C  s2cl 

lw  L  ci +  y  +  y  =  *2  * 

8X  3X  8y 

(4.11) 

with 

X  =  x+/L+  =  x/L  , 

(4.12) 

*  +  1/3 ,  +1/3 

y  =  y  Sc  /L 

(4.13) 

*  mv  1/3 

“  =  *2  Sc 
u 

(4.14) 

The  boundary  conditions  are 

c±  (X,  0)  =  0  0  <  X  <  1 

c1  (X,  “)  =  0 
c1  (0,  y*)  =  0  . 


(4.15) 
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The  transfer  coefficient  can  also  be  represented  as  the  sum  of 
an  average  and  time  dependent  term 

K  *  K  +  k  (4.16) 

Again,  because  of  the  linear  assumption, 

k  =  k  elut  ,  (4.17) 


where  k  is  a  complex  amplitude.  For  the  case  where  a  pseudosteady  state 
assumption  can  be  made,  one  obtains  from  (4.7) 


k 

s 


0.807 

3 


2/3 
» _ 

1/3 


(4.18) 


Since 


s  a 


S  C 


1  3c, 


b  J°  3y  v= 


dX  , 


y=0 


it  follows  from  (4.11),  (4.15)  and  (4.18)  that 


J<_  ,  ,  *+2/3. 

£  =  f  (»  L  )  , 


(4.19) 


(4.20) 


where  k  =  k_  +  i  k_  is  complex  and  k  is  real. 

K  1  s 

The  periodic  variation  of  k,  defined  by  (4.17),  is  characterized 
by  an  amplitude, 


^2  ~  2 

(\)  +  (kx)Z 


(4.21) 
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and  by  a  phase  angle 


-1 


k 


I 


0 


tan 


(4.22) 


k, 


where  positive  0  means  k  lags  s  in  phase.  The  term  (4.22)  represents 
the  phase  correction  and 


k 


s 


(4.23) 


A  = 


the  amplitude  correction  that  must  be  applied  to  the  signal  coming  from 
the  probe,  if  the  pseudosteady  state  relation  (4.7)  is  used  to  calculate 
S(t)  from  the  time  varying  signal. 

Values  of  A  and  0  were  obtained  by  numerically  solving  (4.11)  with 

a  scheme  similar  to  that  used  by  Fortuna  &  Hanratty  (1971). 

2  *  +2/3 

Plots  of  1/A  and  phase  lag,  0,  are  plotted  versus  to  L  in 

Figures  4.2  and  4.3.  The  much  stronger  effect  of  the  frequency  response 

of  the  scalar  boundary- layer  on  the  phase  lage,  0,  is  immediately  evident. 

^  i  2  /  q  2 

At  co  L  '  =1  the  amplitude  correction  1/A  is  only  about  1.07;  however, 
at  this  same  condition  the  phase  lag  correction  will  be  close  to  16°. 

The  principal  unconfirmed  assumption  in  the  above  analysis  is  that 
the  fluctuating  scalar  field  is  described  by  a  linearized  form  of  the 
conservation  equation.  In  order  to  verify  the  numerical  results  obtained 
in  Figures  4.2  and  4.3,  experimental  tests  were  carried  out  for  four 
different  sizes  of  mass  transfer  probes. 


Table  4.1  gives  the  dimensions  of  the  electrodes  and  compares  the 
estimated  thickness  of  the  concentration  boundary-layer  with  the  thickness 
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Figure  4.2  Correction  factor  of  amplitude  for  frequency  response 
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Figure  4.3  Phase  lag  due  to  frequency  response 
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Table  4.1 

Dimensions  of  the  electrodes  and  experimental  range 
for  the  test  of  frequency  response 


Electrode 

L  (cm) 

L+ 

<5+ 

c 

<5+ 

s 

*  +2/3 
oj  L 

1 

.00762 

.352-. 634 

.22-. 27 

2. 5-7. 6 

.23-1.4 

2 

.0127 

.587-1.06 

.27-. 32 

2. 5-7. 6 

.32-1.9 

3 

.0254 

1.17-2.11 

.33-. 41 

2. 5-7. 6 

.51-3.1 

4 

.0406 

1.88-3.38 

.39-. 48 

2. 5-7. 6 

.69-4.2 

Re- =  15,000  -  29,400 
f  =  .325  -  .900  1/sec 


W  =  . 102  cm 


of  the  Stokes  layer  in  the  range  of  variables  covered  by  the  experiments. 
Here  L  is  the  length  of  the  electrode  in  the  direction  of  mean  flow  and 
W,  the  width  of  the  electrode. 

The  thickness  of  the  concentration  boundary-layer  is  defined  as 

,+1/3 

6Jx=L  =  2,92  1/3  *  ’  (4-24) 

Sc  u 


and  the  Stokes  layer  thickness  is  described  in  (3.32). 

«  "b  •  “b 

As  shown  in  Table  1,  6  is  one  order  of  magnitude  less  than  6 

^  s 

and  the  viscous  sublayer  thickness  of  6^  =  5.  Consequently  the 
assumption  of  a  linear  variation  of  the  velocity  within  the  concentra¬ 
tion  boundary  layer  is  valid. 

The  analog  signals  from  the  four  probes  were  sampled  simultaneously 
and  were  processed  in  a  manner  described  in  section  (4.5).  Figures  4.4 
and  4.6  show  the  uncorrected  amplitudes  and  phases  characterizing  the 
variation  of  the  wall  shear  stress  measured  with  the  different  probes. 
Here  the  amplitude  |x  ]  is  plotted  as 


t(o) 


(4.25) 


so  as  to  normalize  effects  of  changes  in  the  amplitude  of  the  induced 
flow  oscillations.  At  low  values  of  the  dimensionless  frequencies, 
there  is  reasonably  close  agreement  among  the  uncorrected  measurements 
from  the  different  probes,  indicating  the  validity  of  the  pseudosteady 
state  assumption.  However,  at  high  large  differences  are  noted. 

Figures  4.5  and  4.7  show  the  results  obtained  by  correcting  the 
measurements  in  Figures  4.4  and  4.6,  using  the  factors  in  Figures  4.2 
and  4.3.  The  good  agreement  now  shown  among  the  measurements  from  all 
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Figure  4.4  Amplitude  of  wall  shear  stress  measured  by  different 
probes  without  a  correction  for  frequency  response 
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Figure  4.5  Amplitude  of  wall  shear  stress  measured  by 
different  probes  with  a  correction  for 
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Figure  4.6  Phase  lag  of  wall  shear  stress  (relative  to  the  trigger  pulse) 
measured  by  different  probes  without  a  correction  for  frequency 
response 
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Figure  4.7  Phase  lag  of  wall  shear  stress  (relative  to  the  trigger  pulse)  measured 
by  different  probes  with  a  correction  for  frequency  response 
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four  probes  provides  support  for  the  method  of  correcting  for  frequency 
response  presented  in  the  numerical  analysis. 

The  effects  of  the  amplitudes  of  flow  oscillation  on  the  amplitude 
of  shear  stress  at  the  wall  are  shown  in  Figure  4.8.  Here  a  is  the  ratio 
of  the  amplitude  of  velocity  variation  to  the  mean  velocity  at  the  pipe 
center.  The  amplitude  of  shear  stress  at  the  wall  is  found  to  increase 
linearly  with  the  increase  of  the  amplitude  of  the  imposed  flow 
oscillation.  If  the  amplitude  | x w |  /  is  plotted  as  t  (0),  the  same 

results  are  obtained  for  different  amplitudes  of  the  imposed  oscillations. 
This  supports  the  linearization  assumption  used  in  deriving  (4.11). 

A  number  of  restrictions  must  be  met  if  wall  transfer  probes  are 
to  be  used  properly  in  a  turbulent  flow  with  imposed  oscillations. 

These  are  summarized  below. 

The  concentration  boundary-layer  should  be  thinner  than  the  viscous 
sublayer  if  (4.2)  is  to  be  satisfied.  This  requires  that 


+  t+1/3 

6cIx+=L+  =  2,92  ~T/1 

S  c 


5 


or  that 


L+  <  5 


Sc 


(4.26) 


The  presence  of  imposed  oscillations  adds  the  additional  restriction 
that  6^  be  less  than  the  thickness  of  the  Stokes  layer: 


+  t+  < 
x  =L 


or 
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Figure  4.8  Effect  of  the  amplitude  of  flow  oscillation  on 
the  amplitude  of  wall  shear  stress 
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-f  +  -3/9 

L  <  0.11  (a>  )  1  Sc.  (4.27) 

The  neglect  of  the  effects  of  diffusion  in  the  flow  direction  (Hanratty 
&  Campbell,  1983)  requires  that 

+  -1/2 

L  >  14  Sc  '  (4.28) 

if  errors  are  to  be  less  than  5  percent. 

Restrictions  (4.26),  (4.27)  and  (4.28)  must  be  satisfied  if  the 
analysis  of  the  time  response  of  the  scalar  boundary-layer  presented 
in  this  work  is  to  be  used.  The  restriction  imposed  by  the  linearization 
assumption  made  in  the  analysis  is  not  as  well  defined.  The  experiments 
described  in  this  study  had  ratios  of  the  amplitude  of  oscillation  in 
the  wall  shear  stress  to  the  time  averaged  wall  shear  stress,  |x  |  /  , 

as  high  as  0.6.  The  successful  experience  in  correcting  these  measure¬ 
ments  suggests  that  (4.28)  might  be  more  restrictive  than  the  linearization 
assumption,  regarding  the  size  of  oscillations  that  can  be  measured.  If 
the  oscillations  become  so  large  that  during  portions  of  the  cycle, 
convection  is  no  longer  completely  dominant  over  streamwise  diffusion 
errors  can  arise.  From  (4.28) 


1  -■ 


w ' 


1/2 


>  14  Sc 


-1/2 


w 


or 


\t  1  ] 

14 

1  w 

T 

V.  w  J 

*> 

Sc1/2L+ 

V.  / 

(4.29) 


If  restrictions  (4.26-4.29)  are  met.  Figures  4.2  and  4.3  can  be 


used  to  correct  for  the  frequency  response  of  the  concentration 
boundary-layer . 

The  numerical  analysis  suggests  that 


+Sc1/3L+2/3 


^  1 


(4.30) 


if  it  is  desired  that  the  correction  on  the  amplitude  is  less  than  5%. 
This  implies  that  it  might  be  necessary  to  correct  measurements  of 
turbulent  spectra  in  the  high  frequency  range.  Fortuna  &  Hanratty  (1971) 
did  this  and  found  that  the  intensity  of  turbulent  shear  stress  fluctu— 
ations  at  the  wall  is  0.349  rather  than  the  value  of  0.302  obtained  by 
using  pseudo-steady  state  relation  (4.7).  The  criterion  for  avoiding 
corrections  of  the  phase  angle  is  more  restrictive: 


ui+Sc1/3l+2/3  <  0>3 


(4.31) 


for  an  error  of  less  than  5°. 

The  strong  effect  of  Schmidt  numbers  is  evident.  Consider  for 
example  Sc  =  1000,  which  might  be  realized  with  an  electrochemical  probe. 
Restrictions  (4.26)  and  (4.28)  indicate  0.44  <  L+  <  5000. 

If  L  =  0.44  is  used  in  (4.27)  it  is  found  that  the  thinness  of 
the  Stokes  layer  is  not  a  limitation  in  that  frequenceis  as  high  at 
w+  =  40  can  be  studied  without  violating  (4.27).  From  (4.31)  it  is 
found,  with  L+  =  0.44,  that  frequencies  as  high  as  w+  =  0.052  can  be 
used  without  correcting  for  the  frequency  response  of  the  scalar  boundary- 
layer.  For  reference,  this  is  approximately  equal  to  the  median  frequency 
of  the  turbulent  velocity  fluctuations  close  to  a  wall  or  to  the  bursting 
frequency.  This  means  that  the  phase  corrections  for  the  frequency 


67 

response  of  the  scalar  boundary-layer  will  have  to  be  taken  into  account 
in  experiments  with  imposed  oscillations,  unless  the  oscillation  has  a 
much  lower  frequency  than  the  turbulent  velocity  fluctuations. 

4.4  Fabrication  of  Rectangular  Electrode 

A  rectangular  electrode  was  built  in  a  plexiglass  block  which  was 
then  inserted  into  the  wall  of  the  test  section.  The  procedure  of 
manufacture  is  shown  in  Figure  4.9.  A  block  of  plexiglass  is  machined 
to  the  configuration  as  shown  in  Figure  4.9(a).  The  depth  of  the  groove 
at  one  end  is  the  thickness  of  the  platinum  sheet,  which  is  the  desired 
length  of  the  electrode.  At  the  other  end  the  depth  of  the  groove  is 
greater,  so  that  it  can  accommodate  the  junction  of  the  platinum  sheet 
and  a  copper  wire.  The  platinum  sheet  is  glued  in  the  groove  with 
epoxy  (Techkits  A-12).  Another  piece  of  plexiglass  is  glued  on  the  top 
of  the  sheet  as  shown  in  Figure  4.9(b).  After  curing  of  the  epoxy,  the 
block  is  machined  to  the  configuration  shown  in  Figure  4.9(c)  and 
inserted  into  a  slit  which  is  cut  parallel  to  the  pipe  axis  in  the  wall 
of  the  test  section.  The  block  is  glued  with  epoxy  in  the  slit  and 
allowed  to  protrude  a  little  to  ensure  the  electrode  is  fully  exposed 
and  yet  not  leave  an  inordinate  amount  of  the  electrode  to  be  abraded 
off.  The  surface  of  the  electrode  is  sanded  and  polished  progressively 
with  wetted  400A,  600A  and  4/0  emery  paper.  Finally  polishing  is  done 
with  Mirror  Glaze  Plastic  Cleaner  and  Polish.  In  this  manner,  the 
rectangular  electrode  can  be  orientated  with  its  longer  side  perpendic— 
ulst  to  the  axial  flow  direction  and  mounted  flush  with  the  wall. 

4 . 5  Differential  Pressure  Measurement 

Since  there  was  a  slight  difference  between  the  phase  and  the 
amplitude  of  the  flow  oscillation  in  the  test  section  and  the  motion 
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Figure  4.9  Fabrication  of  rectangular  electrode 


of  the  pulse  generator,  it  was  necessary  to  measure  directly  the  flow 
oscillation  in  the  test  section.  This  was  done  with  a  differential 
pressure  transmitter  (Viatran,  Model  704-115)  having  a  full  scale 
±  2.54  cm  water  column  corresponding  to  ±  5  volts  output.  Two  pressure 
taps,  437.4  cm  apart,  were  located  on  the  two  ends  of  the  test  section. 

Fora  fully  developed  pipe  flow,  p  is  dependent  of  x  and  t,  so  the 
pressure  gradient  is  only  a  function  of  time.  That  is. 


9p  _  Ap 
3x  Ax 


(t) 


(4.32) 


This  measured  time  variation  of  the  pressure  gradient  was  used,  along 
with  equation  (1.4),  to  calculate  the  time  variation  of  the  velocity 
at  the  center  of  the  pipe. 

The  dynamic  response  of  the  pressure  transmitter  was  determined 
in  the  following  way:  The  flow  loop  was  first  filled  with  the  solution. 
Then,  the  flow  in  the  19.4  cm  pipe  was  oscillated  by  the  motion  of  the 
piston  with  diaphragm  valve  3  completely  closed.  Because  of  the 
incompressibility  of  the  fluid,  the  motion  in  the  test  section  was  in 
phase  with  the  piston  and  the  amplitude  of  the  pressure  gradient 
variation  could  be  calculated  from  the  piston  stroke  length  and  frequency. 
The  instanteous  differential  pressure  signals  were  sampled  and  processed 
as  described  in  the  next  section.  Several  frequencies  and  stroke  lengths 
were  tested.  It  was  found  that  the  phase  lags  of  the  pressure  trans¬ 
mitter  was  dependent  on  frequency  as  shown  in  Figure  4.10.  It  was  9° 

for  a  frequency  of  0.325  Hz  and  25°  for  0.625  Hz.  The  amplitude 

2 

correction  was  dependent  on  L  f  ,  where  L  is  the  stroke  length  of  the 

s  s 

pump.  Figure  4.11  compares  the  measured  amplitude  with  that  calculated 
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Figure  4.10  Phase  correction  for  pressure  transmitter 
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Figure  4.11  Amplitude  calibration  of  pressure  transmitter 
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were  taken  for  each  run.  A  sampling  program  was  written  for  this 
purpose  and  is  given  in  Appendix  III. 

In  the  measurement  of  pulsating  flow,  an  optointerruptor  on  the 
pulse  generator  was  used  to  provide  a  pulse  of  4  volts  when  the  piston 
reached  a  certain  position  in  each  cycle.  This  pulse  served  as  a 
phase  reference.  The  time  between  the  two  consecutive  pulses  is  the 
period  of  the  oscillation.  It  was  measured  with  a  programmable  timer 
(Berkeley,  6401)  and  sampled  through  a  real  time  clock  board  in  the 
computer.  Drift  of  the  period  of  different  cycles  was  found  to  be  less 
than  0.1%.  For  each  cycle  of  oscillation,  this  pulse  initiated  data 
sampling.  Thirty  two  or  sixty  four  samples  per  cycle  were  taken  in 
equal  time  intervals  and  a  total  of  500  periods  of  data  were  collected 
on  the  disc  of  the  computer  in  each  experiment.  At  each  sampling 
point  the  probes  and  the  pressure  transmitter  were  sampled  at  a  high 
enough  speed  so  that  the  measurements  for  all  practical  purposes , 
could  be  considered  simultaneous.  A  program  which  controls  this  sampling 
process  is  given  in  Appendix  III. 

The  digitized  data  were  recorded  on  the  disc  in  a  binary  format 
in  order  to  save  storage  space.  During  data  processing,  the  binary 
data  were  first  converted  to  voltages  in  decimal  format.  These  voltage 
data  correspond  to  the  output  voltages  of  the  electrode  circuit  and  the 
pressure  transmitter  at  each  sampling  instant.  The  output  voltage  of 
the  electrode  circuit  is  related  to  the  velocity  gradient  at  the  wall: 
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B(V 


out 


-  v  ) 

app 


(4.33) 


according  to  (4.1)  and  (4.7),  where  VQut  is  the  output  voltage  and  V  , 
the  applied  voltage  in  the  electrode  circuit.  Parameter  3  is  a  constant 
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for  each  run.  Therefore,  the  periodic  variation  of  wall  shear  stress 
normalized  by  the  time  averaged  value  was  obtained  as 
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The  intensity  of  the  turbulent  velocity  gradient  fluctuations  at 
the  wall  was  given  as 
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(4.35) 


and  the  intensity  of  the  phase  averaged  turbulent  fluctuations,  as 
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The  output  voltage  of  the  pressure  transmitter  is  proportional  to 
the  instantaneous  pressure  drop  (2.54  cm  water  column  per  volt).  The 
time  averaged  pressure  drop  was  calculated  by  the  Blasius  formula  rather 
than  from  the  measured  values,  because  serious  errors  were  caused  by 
the  large  zero  drift  of  the  transmitter  and  the  small  value  of  time  mean 
pressure  drop. 

Phase  averages  and  time  averages  were  calculated  on  the  computer 
in  the  following  manner: 
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where  <  E  >  are  phase  averaged  data  in  one  period  and  n,  the  number  of 
periods  to  be  averaged.  Time  averaging  was  done  by  averaging  the  phase 
averaged  data  over  one  period: 
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where  m  is  the  number  of  samples  in  one  period.  The  oscillating 
component  was  then  obtained  as 

E  =  <  E  >  -  E  .  (4.39) 

Only  the  first  harmonic  of  the  oscillation  was  considered,  so  the 
amplitude  and  the  phase  of  the  velocity  gradient  at  the  wall  were 
obtained  by  a  least-square  fit  of  the  phase  averaged  data  with  a  cosine 
function. 

The  spectrum  of  the  turbulent  fluctuations  of  the  velocity  gradient 
at  the  wall  was  obtained  by  a  Fast  Fourier  Transformation  (Bendat  and 
Piersol  1971).  The  analog  singals  from  the  wall  probe  were  filtered 
by  an  analog  low  pass  filter  with  a  cut-off  frequency  20  Hz  before  they 
were  sampled  by  the  computer  at  a  rate  of  40  Hz.  The  cut-off  frequency 
was  much  higher  than  the  medium  frequency,  f+  =  .01.  Therefore, 
turbulent  fluctuations  with  frequencies  higher  than  the  cut-off  frequency 
make  a  negligible  contribution  to  the  total  energy  (Hanratty,  et_  al. , 
1977).  The  actual  sample  size  for  each  of  FFT  was  512.  A  total  number 
of  20,480  of  data  points  were  used,  i.e.,  40  realizations  were  averaged 
to  give  an  estimated  spectrum.  The  spectrum  was  normalized  so  that  the 


area  under  the  curve  is  unity. 
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4. 7  Experimental  Procedure  and  Range 

Before  doing  an  experiment  the  solution  was  circulated  between  the 
filter  and  the  lower  storage  tank  for  eight  hours  in  order  to  remove 
small  particles  or  dirt.  Then  the  solution  was  pumped  into  the  upper 
storage  tank  and  the  loop  was  filled  through  the  downflow  pipe  line. 
However  the  piping  above  the  solution  level  in  the  tank  couldn't  be  filled 
with  the  solution.  When  the  centrifugal  pump  started,  the  air  inside 
this  portion  of  the  piping  was  pushed  into  the  tank  and  part  of  it  was 
entrained  into  the  flow  loop.  Another  origin  of  the  air  bubbles  in  the 
fluid  was  the  decrease  of  solubility  of  the  air  with  the  increase  of 
solution  temperature  when  the  solution  was  circulated  through  the  loop. 
These  bubbles,  sometimes,  stuck  on  the  wall  of  the  test  section  and 
made  measurement  impossible.  It  took  a  long  time  to  get  rid  of  the 
bubbles.  The  flow  was  turned  on  at  a  high  flow  rate  and  then  turned 
off  to  settle  for  an  hour.  This  process  was  repeated  several  times. 
Sometimes,  it  was  necessary  to  run  the  flow  for  a  long  time  to  increase 
the  temperature  of  the  solution  above  the  experimental  temperature. 

Then  the  flow  was  allowed  to  settle  for  an  hour.  Afterwards,  the  loop 
was  run  again  and  the  flow  temperature  was  cooled  down  to  the  experimental 
temperature  25  ±  0.1  °C  by  adjusting  the  cooling  water  through  the  cooling 
coil.  This  was  a  very  tedious  process,  but  it  was  very  important  for 
obtaining  satisfactory  data. 

Mean  flow  rate  was  adjusted  by  diaphragm  valve  3  and  was  measured 
by  the  vortex  shedding  flow  transmitter.  For  steady  flow  measurements, 
at  least  thirty  minutes  for  stabilizing  the  flow  was  necessary  before 
taking  data.  For  the  measurements  of  the  pulsating  flow,  the  frequency 


of  the  oscillation  was  adjusted  to  the  desired  value  and  was  kept 
constant  for  at  least  500  periods  before  sampling  the  data  in  order 
to  reach  stationary  oscillation. 

The  polarization  curves  of  the  electrode  in  the  experimental 
range  of  Reynolds  numbers  were  measured,  as  shown  in  Figure  4.14.  These 
suggested  the  choice  of  -0.4  volts  as  the  applied  voltage  to  the  test 
electrode. 

The  wave  form  of  the  oscillation  was  monitored  by  the  differential 
pressure  signal  on  the  oscilloscope.  If  the  pressure  signal  was  distorted 
from  a  sinusoidal  variation,  it  was  likely  to  that  air  was  in  some  portion 
of  the  system.  Data  were  not  taken  until  the  pressure  signal  showed  a 
sinusoidal  variation.  The  period  of  each  cycle  was  indicated  by  the 
programmable  timer.  If  the  drift  of  the  period  was  greater  than  0.1% 
of  its  value  in  the  experiment,  the  data  were  disregarded. 

The  experiments  were  designed  to  cover  a  range  of  oscillation 
frequencies  around  the  medium  frequency  of  the  turbulent  velocity 
fluctuations  in  the  vicinity  of  the  wall,  i.e.,  u)+  =  2tt  x  0.01.  This 
was  realized  by  fixing  the  frequency  of  the  imposed  oscillation  at 
0.625  Hz  or  0.325  Hz  and  varying  the  Reynolds  number  of  the  mean  flow 
from  15,000  to  70,000.  The  amplitude  of  the  central  velocity  variation 
was  10%  of  the  mean  velocity  for  most  of  the  experiments.  However, 
for  high  Reynolds  numbers,  an  amplitude  of  less  than  10%  had  to  be  used 
because  of  the  limited  capacity  of  the  pulse  generator.  An  oscillation 

m 

amplitude  of  5%  was  also  used  at  Reynolds  numbers  of  15,000  and  20,000 
with  a  frequency  0.625  Hz  in  order  to  examine  the  assumption  of  a 
linear  response  of  the  wall  shear  stress  to  the  imposed  oscillation. 
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Figure  4.14  Polarization  of  electrode  (L  =0.00762  cm) 


The  reason  for  using  a  small  amplitude  was  to  eliminate  nonlinear 
effects  and,  therefore,  to  make  the  analysis  of  the  results  easier. 
The  experimental  conditions  are  listed  in  Table  4.2.  Several 
dimensionless  groups  of  frequency  are  also  calculated  for  comparison 
with  studies  of  other  researchers. 


Table  4.2 

Experimental  Range 
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8.7 
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149 

99 
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0.0074 
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77 

4.4 

0.0045 
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0.325 

11.9 

0.1 

149 

54 

3.0 

0.0023 
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0.325 

12.7 
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149 

40 

2.2 

0.0013 
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0.325 

12.7 
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149 

35 

1.8 

0.0009 

60,000 

0.325 

12.7 

0.052 

149 

29 

1.5 

0.00065 
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0.325 

12.7 

0.044 

149 

26 

1.3 

0.00050 
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0.625 

3.2 

0.1 

207 

190 

11.4 

0.014 

20,000 
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4.3 

0.1 
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8.53 

0.0086 

30,000 
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6.4 

0.1 

207 

104 

5.69 

0.0042 

40,000 

0.625 

8.5 

0.1 

207 

77 

4.27 

0.0026 

50,000 

0.625 

10.7 

0.1 

207 

66 

3.41 

0.0017 

60,000 

0.625 

12.8 

0.1 

207 

57 

2.84 

0.0013 

70,000 

0.625 

12.7 

0.086 

207 

49 

2.44 

0.00095 

CHAPTER  5 .  RESULTS 


5.1  Numerical  Results 

The  numerical  scheme  was  first  tested  by  comparing  the  numerical 
result  for  the  quasilaminar  model  A  with  the  asymptotic  solution  for 
laminar  pulsating  flow  at  high  frequency.  Figure  5.1  shows  velocity 
Profiles  of  the  oscillating  component,  normalized  by  the  amplitude  of 
velocity  at  the  pipe  center,  at  eight  phases  of  a  cycle  at  Q  =  206. 
Figures  5.2  and  5.3  give  the  amplitude  profiles  and  the  phase  profiles 
of  the  oscillating  velocity  component  for  the  two  frequencies  used  in 
the  present  study.  The  solid  lines  represent  the  numerical  results  of 
Model  A  while  the  symbols  are  the  asymptotic  results  obtained  from  the 
analytical  solution  for  high  frequencies.  Agreement  between  the  two 
is  very  good.  This  indicates  that  the  numerical  scheme  used  in  this 
study  is  suitable  for  further  studies  of  the  four  models  described  in 
Chapter  3. 

Figures  5.4  and  5.5  give  the  amplitude  profiles  and  the  phase 
profiles  of  oscillation  velocity  calculated  from  four  different  models 
at  Re  =  15,000  and  f  =  0.625  Hz.  This  corresponds  to  the  high  co+ 
extreme  covered  in  the  experiments.  It  is  seen  that  the  variation  of 
the  oscillating  velocity  is  confined  to  a  thin  layer  close  to  the  wall 
about  0.01  of  the  pipe  radius,  which  corresponds  to  y+  =  5,  i.e.,  the 
thickness  of  viscous  sublayer.  The  four  models  behave  very  similarly 
because  of  the  negligible  effects  of  turbulence  in  the  viscous  sublayer. 
In  the  range  of  0.01  <  Y  <  0.04  the  four  models  have  slight  differences 
because  they  use  amounts  of  oscillation  induced  eddy  viscosity  in  this 
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Figure  5.1  Oscillation  velocity  profiles  at  ft  =  206 
( —  numerical  solution  of  Model  A 
symbol:  asymptotic  solution  of  laminar 

pulsating  flow  at  high  frequency) 
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Figure  5.2  Amplitude  profiles  of  oscillation  velocity 
( —  numerical  solution  of  Model  A 
symbol:  asymptotic  solution  of  laminar 
pulsating  flow) 
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Figure  5.3  Phase  profiles  of  oscillation  velocity 
( —  numerical  solution  of  Model  A 
symbol:  asymptotic  solution  of  laminar 

pulsating  flow) 
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Figure  5.4  Amplitude  profiles  of  oscillation  velocity  for 

different  models  at  Re  =  15,000  and  f  =  0.625  Hz 
((0/15  =  0.014) 


86 


Y*10 

»  — *  ■  ■  *  .  .  .  .  I 

5  10  15  20  25  30  35  40  45 

Y+ 


Figure  5.5  Phase  profiles  of  oscillation  velocity  for 

different  models  at  Re  =  15,000,  f  =  0.625  Hz 
(m+/15  =  0.014) 
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range.  The  models  C  and  D  give  larger  oscillation  induced  eddy 
viscosities  than  models  A  and  B,  so  that  the  over  shoots  of  amplitude 
for  models  C  and  D  are  less  then  models  A  and  B.  The  difference 
between  model  C  and  model  D  is  caused  by  the  pressure  gradient  effect 
considered  in  model  D.  In  the  range  of  Y  >  0.04  the  results  of  the 
four  models  become  closer  again.  This  is  because  inertial  effects  are 
so  large  in  this  region  that  differences  in  eddy  viscosities  predicted 
by  the  four  models  have  negligible  influence  on  the  solution. 

Figures  5.6  and  5.7  show  the  amplitude  profiles  and  the  phase 
profiles  calculated  from  the  four  models  at  Re  =  60,000  and  f  =  .325  Hz. 
The  profiles  of  model  A  are  very  similar  to  that  shown  in  Figures  5.4 
and  5.5,  except  that  the  hump  in  the  amplitude  profile  moves  slightly 
away  from  the  wall  due  to  the  lower  frequency  in  this  case.  However, 
the  results  from  models  B,  C  and  D  are  quite  different.  At  a  Reynolds 
number  60,000  the  viscous  sublayer  is  much  thinner  than  at  Re  =  15,000. 
Although  the  spatial  variation  of  oscillation  velocity  is  in  a  thin 
layer  of  Y  =  0.02,  the  turbulent  effects  on  the  flow  oscillation  are 
no  longer  negligible  because  this  layer  corresponds  to  y+  «  25.  It 
is  seen  in  Figure  5.7  that  the  phase  of  velocity  gradient  at  the  wall 
is  very  sensitive  to  the  different  turbulent  models.  In  the  range  of 
Y  >  0.08  the  inertial  force  controls  the  flow  oscillation  and  the 
different  models  show  similar  results. 

From  the  comparison  of  different  turbulent  models,  it  is  clearly 
seen  that  the  modeling  of  turbulence  in  the  region  of  y+  <  40  is  the 
key  to  understanding  turbulent  pulsating  flow  at  high  oscillation 
frequencies  and  that  the  relative  thickness  of  the  stokes  layer  and  the 
viscous  sublayer  is  a  measure  of  the  importance  of  turbulent  effects  on 
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Figure  5.6  Amplitude  profiles  of  oscillation  velocity  for 
different  models  at  Re  =  60,000,  f  =  0.325  Hz 
(m+/15  =  0.00065) 
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Figure  5.7  Phase  profiles  of  oscillation  velocity  for  different 
models  at  Re  «  60,000,  f  -  0.325  Hz  (o)+/15  =  0.00065) 
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the  flow  oscillation.  It  is  also  found  that  the  wall  shear  stress  is 
very  sensitive  to  the  turbulent  modeling.  These  types  of  numerical 
results  will  be  compared  with  experimental  data  in  section  5.3. 

5.2  Experimental  Results 

(a)  Instantaneous  velocity  gradient  at  the  wall 

Figure  5.8  shows  typical  measurements  of  the  streamwise  components 
of  the  instantaneous  velocity  gradient  at  the  wall  for  a  Reynolds  number 
30,000.  In  Figure  5.8(a),  the  time  averaged  value  is  subtracted  from 
the  instantaneous  value,  so  the  signal  represents  a  summation  of  the 
oscillation  component  and  turbulent  fluctuations.  Since  turbulent 
fluctuations  at  the  wall  are  large,  the  organized  part  of  velocity  gradient 
is  buried  in  random  fluctuations.  However,  the  periodic  variation  can  still 
be  distinguished  when  the  signal  is  compared  to  the  ticks,  representing 
the  time  period  of  the  flow  oscillation,  in  Figure  5.8(a).  It  is 
observed  that  the  period  of  the  organized  variation  in  5.8(a)  is  of  the 
same  order  as  the  period  of  the  characteristic  turbulent  fluctuations  in 
a  steady  flow  at  the  same  Reynolds  number  30,000,  shown  in  Figure  5.8(c). 
This  is  expected  since  the  imposed  oscillation  frequency,  0.625  Hz,  was 
chosen  to  match  the  characteristic  frequency  of  the  turbulent  fluctuations 
in  the  vicinity  of  the  wall  at  Reynolds  number  30,000. 

In  Figure  5.8(b),  phase  averaged  values  are  subtracted  from  the 
instantaneous  signals,  i.e.,  the  oscillation  component  is  removed  from 
5.8(a).  The  signal  in  Figure  5.8(b)  represents  the  turbulent  fluctuations 
of  the  velocity  gradient  at  the  wall  in  the  pulsating  flow.  A  comparison 
5.8(b)  with  5.8(c)  shows  no  major  differences.  This  indicates  that  for 
imposed  oscillations  of  small  amplitude,  the  interaction  between  the 
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Figure  5.8  Instantaneous  signals  from  wall  probe  at  Re  =  30,000 

(a)  with  oscillation,  f  =  0.625  Hz  (w*7l5  =  0.0^2) 

(b)  with  oscillation,  f  =  0,625  Hz,  but  subtracting  oscillation 
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turbulence  and  the  imposed  oscillation  is  not  clearly  seen,  even  though 
their  frequencies  are  similar. 

Power  spectra  of  the  instantaneous  signals  are  shown  in  Figure  5.9. 

The  sharp  peak  at  a  frequency  0.63  Hz  in  Figure  5.9(a)  represents  the 
contribution  of  the  imposed  oscillation  to  the  total  energy.  Figure  5.9(b) 
shows  the  spectrum  of  turbulent  fluctuations  of  velocity  gradient  at  the 
wall  after  removing  the  phase  averaged  values  from  the  instantaneous 
signals.  Figures  5.9(c)  is  the  spectrum  of  turbulent  fluctuations  of 
the  velocity  gradient  at  the  wall  in  steady  flow  at  a  Reynolds  number 
30,000.  No  major  difference  is  observed  between  Figures  5.9(b)  and  4.9(c). 
This  would  seem  to  suggest  that  the  turbulent  properties  are  not  affected 
by  the  imposed  oscillation. 

The  peaks  at  frequencies  13  Hz  and  3  Hz  are  considered  to  be 
caused  by  some  unknown  disturbance  in  the  flow  system,  because  the  peaks 
appear  at  the  same  frequencies  in  the  steady  flow  at  different  Reynolds 
numbers. 

(b)  Time  averaged  values  of  the  velocity  gradient  at  the  wall 

Figure  5.10  shows  the  measured  time  averaged  values  of  the  velocity 
gradient  at  the  wall  at  different  flow  conditions.  It  is  observed  that 
the  imposed  oscillation  has  no  effect  on  the  time  mean  velocity  gradient 
at  the  wall. 

Intensities  of  turbulent  fluctuations  of  the  velcotiy  gradient  at 
the  wall  are  shown  in  Figure  5.11.  The  intensities  are  the  results  of 
long  time  averaging,  and  are  noticed  to  be  independent  of  Reynolds 
number  and  of  imposed  flow  oscillation.  The  value,  0.35,  agrees  with 
previous  measurements  by  other  researchers  (Sirkar  1969,  Fortuna  and 
Hanratty  1971,  Eckelman  1971). 


with  oscillation  f  =  0.62  Hz  subtract  oscillation  component  steady  state 
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Figure  5.9  Spectrum  of  fluctuations  of  velocity  gradient 
at  the  wall  (Re  =  30,000) 
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Figure  5.10  Time  average  velocity  gradient  at 

the  wall  in  different  flow  conditions 
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Figure  5.11  Time  averaged  intensities  of  turbulent  fluctuations 
of  velocity  gradient  at  the  wall  in  different  flow 
conditions 


The  results  that  have  been  presented  here  indicate  that  time  averaged 
properties  in  a  pulsating  flow  are  not  different  from  those  observed  in  a 
flow  without  pulsations.  These  agree  with  the  findings  by  most  of 
previous  investigators  and  serve  as  a  base  to  separate  the  oscillation 
flow  field  by  simply  subtracting  the  mean  flow  field,  which  is  the  same 
as  the  steady  flow,  from  the  phase  averaged  flow  field. 

(c)  Phase  averaged  values  of  the  velocity  gradient  at  the  wall  in  a 

pulsating  flow 

Figures  5.12  to  5.26  give  a  set  of  phase  averaged  data  measured  at 
different  flow  conditions  in  the  experimental  range.  In  each  figure, 
the  phase  averaged  favorable  pressure  gradient  normalized  by  its  time 
mean  value  is  given  in  (a) ,  the  phase  averaged  vleocity  gradient  at  the 
wall  or  wall  shear  stress  normalized  by  its  time  mean  value  is  shown  in 
(b)  and  the  phase  averaged  intensity  of  turbulent  fluctuations  of  velocity 
gradient  at  the  wall  is  plotted  in  (c) .  The  abscissa  is  the  phase  in  a 
cycle,  referred  to  the  trigger  signal  from  the  piston  as  described  in 
Chapter  4.  The  data,  which  are  the  phase  averaged  results  for  the 
instantaneous  signal  over  500  periods,  have  been  corrected  for  errors 
caused  by  the  frequency  response  of  the  sensors.  In  the  following,  two 
extremes  in  the  experimental  range  will  be  discussed  in  detail. 

Figure  5.12(a)  shows  the  phase  averaged  pressure  gradient  measured 
at  a  Reynolds  number  of  15,400  using  an  imposed  0.625  Hz  oscillation 
which  has  an  amplitude  10%  of  the  central  mean  velocity.  It  is  seen 
that  the  amplitude  of  the  pressure  gradient  oscillation  is  very  large 
in  this  case,  about  eighty  times  its  mean  value.  This  is  caused  by  the 
relatively  high  frequency  of  the  oscillation  and  the  low  time  mean 
velocity.  This  case  represents  one  extreme  in  the  experimental  range, 
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i.e.,  a  low  Reynolds  number  and  a  high  frequency  oscillation.  The  solid 
line  is  the  least  square  fit  of  the  data  with  a  cosine  curve.  It  is 
seen  that  most  of  the  energy  of  the  flow  oscillation  is  contained  in  a 
single  harmonic. 

Figure  5.12(b)  shows  the  phase  averaged  vleocity  gradient  at  the 
wall,  normalized  with  its  mean  value.  The  50%  variation  of  velocity 
gradient  at  the  wall,  compared  with  the  10%  variation  of  the  central 
velocity,  arises  because  the  spatial  variation  of  the  oscillation 
velocity  is  confined  to  a  thin  layer  close  to  the  wall  for  this  high 
frequency  case.  The  data  are  fit  quite  well  with  a  cosine  curve  which 
is  shifted  about  54°  relative  to  the  curve  describing  the  pressure 
gradient.  This  indicates  a  linear  response  of  the  phase  averaged  flow 
field  to  the  imposed  pressure  variation.  A  linear  response  is  also 
indicated  from  experimental  results  shown  in  Figure  5.13,  which  were 
obtained  with  a  5%  variation  of  the  central  velocity.  The  amplitude 
of  pressure  gradient  is  about  fourty  times  that  of  the  mean  and  the 
amplitude  of  velocity  gradient  is  one  half  of  that  indicated  in  Figure 
5.12(b).  The  phase  shifts  of  the  velocity  gradients  are  the  same  for 
both  cases. 

It  is  interesting  that  phase  averaged  intensities  of  turbulent 
fluctuations  are  highly  nonlinear,  as  shown  in  Figure  5.12(c),  even 
though  their  time  average  is  the  same  as  in  steady  flow.  This  nonlinear 
effect  seems  to  be  related  to  the  amplitude  of  the  imposed  oscillation. 
For  example,  the  nonlinear  effect  is  not  as  obvious  in  Figure  5.13(c) 
as  it  is  in  Figure  5.12(c).  It  is  of  interest  to  note  that  the  maximum 
value  of  intensity  occurs  close  to  when  the  pressure  gradient  has  its 
maximum  favorable  value,  and  the  minimum,  close  to  when  the  maximum  of 
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unfavorable  pressure  gradient  occurs.  This  contrasts  with  the  finding 
for  boundary  layer  flows  that  turbulence  near  the  wall  is  suppressed 
by  a  favorable  pressure  gradient  and  enhanced  by  an  unfavorable  pressure 
gradient.  These  results  indicate  that  the  turbulence  does  not  respond 
immediately  to  a  high  frequency  oscillation  of  the  pressure  gradient. 

Figure  5.25  shows  another  set  of  data  at  a  Reynolds  number  60,000 
and  a  frequency  0.325  Hz.  The  amplitude  of  the  central  velocity 
oscillation  is  about  5%.  This  case  represents  another  extreme  in  the 
experimental  range;  that  is,  a  high  Reynolds  number  and  a  low  frequency. 
In  Figure  5.25(a),  the  amplitude  of  the  pressure  gradient  is  seven  times 
its  time  mean  value  and  there  is  some  scatter  of  the  data.  However, 
the  data  still  clearly  show  a  sinusoidal  variation  if  only  the  first 
harmonic  is  considered.  Phase  averaged  data  of  the  velocity  gradient 
at  the  wall  are  shown  in  Figure  5.25(b).  For  the  same  5%  amplitude  of 
central  velocity  variation,  the  responses  at  the  wall  for  the  two 
extremes  of  the  experiments  are  quite  different,  as  can  be  seen  by 
comparing  Figure  5.13(b)  and  5.25(b).  For  example,  the  amplitude  of 
the  velocity  gradient  is  about  6%  and  the  phase  shift  is  67°  in  this 
case  rather  than  about  25%  and  54°  in  the  case  of  low  Reynolds  number 
and  high  frequency.  The  intensity  data  shown  in  5.25(c)  are  more 
scattered.  However,  a  sinusoidal  variation  can  still  be  distinguished. 
It  is  interesting  to  note  that  in  this  case  the  minimum  value  of  the 
intensity  occurs  close  to  the  maximum  favorable  pressure  gradient  and 
the  maximum  occurs  near  the  maximum  of  unfavorable  pressure  gradient. 
This  indicates  that  high  frequency  turbulence  that  occurs  at  large 
Reynolds  number  can  adjust  to  a  slow  oscillation  of  the  pressure 
gradient. 
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Some  general  trends  can  be  seen  from  this  set  of  figures.  For 
a  fixed  oscillation  frequncy  and  amplitude  of  the  central  velocity,  the 
normalized  amplitude  of  pressure  gradient  decreases  as  the  Reynolds 
number  increases;  so  does  the  amplitude  of  wall  shear  stress.  This  is 
very  different  from  the  results  of  laminar  pulsating  flow,  which  depends 
only  on  the  dimensionless  frequency  ft.  This  difference  can  be  attributed 
to  the  strong  effects  of  turbulence  on  the  flow  oscillation.  For  the 
same  Reynolds  number  and  amplitude  of  flow  oscillation,  a  higher 
oscillation  frequency  gives  higher  an  amplitude  of  the  variation  of  the 
wall  shear  stress.  The  amplitude  effect  is  linear  when  it  is  small,  as 
shown  in  Figures  5.12,  5.13  and  4.8. 

At  large  Reynolds  number  the  imposed  oscillation  is  found  not  to  be 
perfectly  sinusoidal.  This  is  attributed  to  the  use  of  a  large  stroke 
length,  which  is  necessary  to  maintain  the  same  amplitude,  as  shown  in 
Appendix  IV. 

It  is  difficult  to  draw  any  quantitative  conclusion  from  the 
scattered  intensity  data.  However,  it  seems  that  certain  qualitative 
trends  can  be  recognized.  For  both  imposed  frequencies  studied,  the 
nonlinear  response  of  the  phase  averaged  turbulent  intensity  of  the 
wall  shear  stress  exists  at  low  Reynolds  number  or  high  <d+.  This  non¬ 
linearity  becomes  smaller  as  Reynolds  number  increases  and  as  the 
amplitude  of  imposed  oscillation  decreases.  At  high  m+  the  variation 
of  the  phase  averaged  intensity  is  roughly  in  phase  with  the  favorable 
pressure  gradient.  At  about  w+/15  =  0.002,  the  periodic  variation  of 
the  intensity  is  hardly  seen,  and  a  further  decrease  of  co"*"  the  variation 
of  the  intensity  in  a  cycle  become  clear  again  but  its  phase  is  roughly 
180°  out  of  phase  with  the  favorable  pressure  gradient. 
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Figure  5.12  Phase  averaged  values  at  Re  =  15,400, 

f  =  0.625  Hz  and  a  =  0.1  (to+/15  =  0.014) 
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Figure  5.13  Phase  averaged  values  at  Rc  p  15,400,  f  =  0.625  Hz 
and  a  =  0.05  (u>+/15  =  0.014) 


INTENSITY 


102 


DEGREE 


Figure  5.14  Phase  averaged  values  at  Re  =  19,900,  f  =  0.625  Hz 
and  a  =  0.1  (w+/15  =  0.0086) 
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and  a  =  0.1  (u)"*"  =  0.0042) 
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and  a  =  0.1  (w+/15  =  0.0025) 
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Figure  5.17  Phase  averaged  values  at  Re  =  49,800,  f  =  0.625  Hz 
and  a  =  0.1  (w+/15  =  0.00174) 
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Figure  5.18  Phase  averaged  values  at  Re  =  57,500,  f  =  0.625  Hz 
and  a  =  0.1  (co+/15  =  0.00135) 
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Figure  5.19  Phase  averaged  values  at  Re  =  70,000,  f  =  0.625  Hz 
and  a  -  0.086  (uj+/15  =  0.00095) 
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and  a  =  0.1  (cj+/15  =  0.0074) 
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and  a  =  0.1  (co+/15  =  0.0045) 
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and  a  =  0.1  (w+/15  =  0.0022) 
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Figure  5.23  Phase  averaged  values  at  Re  =  40,400,  f  =  0.325  Hz 
and  a  =  0.077  (w+/15  =  0.0013) 
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Figure  5.24  Phase  averaged  values  at  Re  =  49,800,  f  =  0.325  Hz 
and  a  =  0.062  (oj+/15  =  0.00091) 
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Figure  5.25  Phase  averaged  values  at  Re  n  605Q0Q,  jj  =  0.325  Hz 
and  a  =  0.05  (w+/15  =  0.00065) 
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Figure  5.26  Phase  averaged  values  at  Re  =  69,800,  f  =  0.325  Hz 
and  a  =  0.044  (oj+/15  =  0.00050) 
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In  summary,  the  phase  averaged  data  indicate  that  the  velocity 
gradient  at  the  wall  responds  to  the  imposed  flow  oscillations  differently 
at  different  flow  conditions. 

5-3  Comparison  of  the  Experimental  Data  and  the  Numerical  Results 

The  velocity  gradient  at  the  wall  shows  a  linear  response  to  sinu¬ 
soidal  imposed  oscillations  of  small  amplitude.  Its  variation  can  be 
characterized  by  an  amplitude  and  a  phase  shift  relative  to  the  imposed 
oscillation.  Figures  5.27  and  5.28  summarize  the  experimental  results 
obtained  under  different  flow  conditions.  The  phase  shifts  are  plotted 
relative  to  the  variation  of  the  central  velocity,  which  lags  the 
pressure  gradient  by  90°,  as  described  in  equation  (1.6).  The  amplitude 
|SXI  /  normalized  by  |u^|  /  u  to  take  account  of  the  effect  of 

using  different  amplitudes  of  oscillation.  The  correlation  is  similar 
to  what  was  used  by  Abrams  and  Hanratty  (1984)  to  characterize  the 
spatial  variation  of  shear  stress  along  a  wavy  surface.  The  abscissa 

j  ^ 

w  /15  corresponds  to  the  a  =  2ttv/Au  used  by  Abrams  and  Hanratty  (1984). 

The  15  is  the  dimensionless  turbulence  convection  velocity  in  the  viscous 

+ 

wall  region.  For  fixed  frequencies,  an  increase  of  -~r-  corresponds  to  a 
decrease  of  Reynolds  number.  Consequently,  the  data  show  that  the 
phases  and  amplitudes  for  fixed  frequencies  are  changing  with  Reynolds 
number.  It  is  noted  that  the  results  obtained  for  the  two  frequencies 
studied  fall  on  the  same  curves  when  plotted  in  the  manner  shown  in 
Figures  5.27  and  5.28.  For  reference,  the  data  obtained  by  Ramaprian 
(1983)  are  also  plotted  in  Figures  5.27  and  5.28.  The  experimental  data 
for  all  runs  are  tabulated  in  Appendix  I. 


116 


o  a: 


o 


m 


(daj&8Q)  pDd~|  asoMd 


Figure  5*27  Phase  lead  (relative  to  the  central  velocity) 
of  velocity  gradient  at  the  wall 
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Figure  5.28  Amplitude  of  velocity  gradient  at  the  wall 


118 


In  Figures  5.27  and  5.28  numerical  results  obtained  with  the 

different  models  described  in  Chapter  3  are  compared  with  the  data.  When 

+ 

■JTj  >  0.01,  the  results  of  all  models  approach  that  predicted  by  a  quasi- 

laminar  model.  Only  Model  D  with  =  -25,  kL  =  3500  predicts  the  rapid 

change  of  phases  and  amplitudes  in  the  range  yr-  -  0.0005  -  0.002.  It  is 

4* 

interesting  to  note  that  this  range  of  is  the  same  range  of  a"*"  in 
which  the  phase  angle  characterizing  the  wall  shear  stress  variation  along 
a  wavy  surface  undergoes  a  rapid  change.  These  results  suggest  that  the 
change  of  phase  in  this  range  for  the  two  different  types  of  flows  can 
be  attributed  to  the  same  relaxation  phenomenon. 

However,  none  of  the  models  are  able  to  predict  the  experimental  data 
in  the  range  of  —  >  0.02.  This  suggests  that  the  models  fail  to  predict 
the  experimental  results  if  the  imposed  oscillation  frequency  is  close 
to  or  greater  than  the  median  frequency  of  the  turbulence  in  the  vicinity 
of  the  wall,  =  0.042.  The  experimental  results  indicate  that  at 
these  high  frequencies  turbulence  is  still  affecting  the  imposed 
oscillations. 

No  reliable  measurements  of  the  oscilation  velocity  profile  were 
obtained  in  this  study.  Attempts  were  made  to  use  a  mass  transfer  probe 
to  measure  the  velocity  in  the  fluid.  Unfortunately,  contamination  of 
the  surface  of  the  sensor  wire  made  reliable  measurement  impossible. 
Moreover,  the  interaction  between  the  wall  and  the  sensor  wire  might 

cause  large  errors  in  the  measurement  close  to  the  wall.  Therefore,  the 

ability  to  predict  the  oscillation  velocity  profile  for  the  four  models 
is  checked  by  comparing  the  numerical  results  with  the  measured  LDV  data 
by  Ramaprian  and  Tu.  The  data  were  taken  at  Re  =  50,000  and  f  =  3.6  Hz 

in  a  water  pipe  flow  with  a  pipe  diameter  of  5  cm.  This  flow  condition 


corresponds  to  the  lowest  a)  in  the  present  experimental  range. 

Figures  5.29  and  5.30  show  these  comparisons  for  the  amplitude  and  the 
phase  of  the  oscillation  velocity  in  the  wall  region.  It  is  seen  that 
Model  D  does  the  best  job  and  agrees  with  the  data  very  well. 
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figure  5.29  Comparison  of  amplitude  profile  among  numerical 
results  and  Ramaprian's  data  (Re  =  50,000, 
f  =  3.6  Hz,  d  =  5  cm) 
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Figure  5.30  Comparison  of  phase  profile  among  the 
numerical  results  and  Ramaprian's  data 
(Re  =  50,000,  f  =  3.6  Hz,  d  =  5  cm) 
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CHAPTER  6.  DISCUSSION 

This  thesis  reports  on  measurements  of  the  time  variation  of  the 
phase  averaged  wall  shear  stress  caused  by  the  imposition  of  controlled 
sinusoidal  oscillations  on  the  flow  of  a  turbulent  fluid  through  a  pipe. 
Results  are  presented  for  a  range  of  dimensionless  frequencies,  uf1",  higher 
than  had  been  presented  by  previous  investigators.  By  using  a  small 
amplitude  of  imposed  oscillation  a  linear  response  is  obtained  whereby 
the  wall  shear  stress  variation  is  described  by  a  single  harmonic  with 
the  same  frequency  as  the  imposed  oscillation  and  with  an  amplitude  which 
varies  linearly  with  the  amplitude  of  the  imposed  oscillation.  The 
measurements  can  therefore  be  characterized  by  a  relative  amplitude  and 
a  relative  phase. 

It  is  found  that  these  measurements,  with  two  different  frequencies 

and  a  range  of  flow  rates,  are  correlated  with  oj+.  This  contrasts  to 

the  characterization  of  turbulent  oscillating  flow  by  coD/u^,  suggested 

by  Jayaraman,  et  al.  (1982),  and  by  a)D  /  u  ,  by  Ramaprian  and  Tu  (1983). 

The  influence  of  frequency  at  high  frequencies  can,  therefore,  be  studied 

either  by  changing  the  flow  rate  or  by  changing  the  frequency. 

+  v2 

The  dimensionless  frequency  is  defined  as  (jj  =  — r^r  /  —  or  as 

*2  a) 

+  2  *  U 
a)  =  2(6  /6  )  ,  where  6  =  v/u  is  the  length  scale  of  the  viscous  sublayer 


v  s 


v 


and  6s  =  /  2 v/o>  is  the  thickness  of  the  Stokes  layer.  Consequently,  for 
>  0.08  (or  o)+/15  >  0.0053)  the  Stokes  layer  is  thinner  than  the  viscous 
sublayer,  and  it  might  be  expected  that  turbulence  will  have  a  small  effect 
on  the  imposed  oscillations.  However,  it  is  to  be  noted  that  measured 
amplitudes  are  much  lower  for  u)+  >  0.08  than  what  is  predicted  for  a  quasi- 
laminar  behavior.  This  result  was  unexpected  and  is  not  yet  understood. 
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Considerable  work  was  done  (Mao  and  Hanratty,  1984)  to  examine  the 
frequency  response  of  the  wall  probes.  The  consequence  of  this  work  is 
that  the  possible  errors  involved  in  correcting  for  frequency  response 
are  far  less  than  the  difference  from  the  quasi-laminar  solution  shown 
in  Figure  5.28.  The  measurements  of  the  oscillation  induced  variation 
in  the  wall  shear  stress  would  suggest  the  possibility  that  the  result 
could  be  due  to  some  nonlinear  effect.  This  explanation  is  ruled  out 
since  the  same  results  have  been  obtained  with  two  different  oscillation 
amplitudes ,  as  shown  in  Figures  5.12  and  5.13.  Furthermore,  measurements 
of  the  velocity  field  very  close  to  the  wall  by  Jayaraman,  et  al.  (1982) 
and  by  Binder,  et  al.  (1982)  also  seem  to  suggest  lower  oscillation 
amplitudes  at  high  frequencies  that  would  be  calculated  by  the  quasi- 
laminar  model.  Consequently,  we  conclude  that  oscillation  induced  varia¬ 
tions  of  the  turbulence  in  the  viscous  wall  layer,  not  taken  into  account 
by  present  theories  of  turbulence,  occur  in  the  viscous  sublayer  at  high 
u  .  However ,  these  variations  appear  to  have  no  measurable  effect  on 
the  time-average  velocity  (Figure  5.10)  or  on  the  time-averaged  turbulence 
properties . 

It  is  to  be  noted  that  the  medium  frequency  of  the  turbulence  is 
w+/15  ~  0.004.  Thus,  the  effect  discussed  above  appears  to  occur  when 
the  imposed  oscillation  is  equal  to  or  greater  than  the  dominant  frequency 
of  the  turbulence.  This  result  seems  consistent  with  the  suggestion  by 
Mizushina,  et  al.  (1973)  that  significant  changes  in  the  turbulence 
structure  occur  when  the  frequency  of  the  imposed  oscillation  is  the  same 
as  the  bursting  frequency  of  turbulence. 


+ 

For  —  <  0.0021  the  dimensionless  Stokes  layer  thickness  is  greater 
than  eight  wall  units,  5+  -  8.  One  would  not  expect  a  quasi-laminar 
model  to  be  applicable.  The  strightforward  application  of  a  turbulence 
model  developed  for  flow  over  a  flat  surface  (Model  C)  is  not  capable 
of  describing  the  experimental  resutls,  and,  in  particular,  the  sharp 
change  of  the  phase  angle  over  a  rather  narrow  range  of  However, 

the  use  of  a  model  (Model  D)  which  argues  that  flow-induced  oscillations 
in  the  pressure-gradient  can  enhance  (unfavorable  pressure  gradient)  or 
dampen  (favorable  pressure  gradient)  turbulence  in  the  viscous  wall 
region  (y+  <  40)  does  appear  to  describe  the  results.  In  fact,  the  para¬ 
meters,  k^  =  -25  and  k^  =  3500,  chosen  to  fit  the  experimental  results 
are  very  close  to  those  used  by  Loyd,  Moffat  and  Kays  (1970)  to 
describing  boundary-layer  flows  and  by  Abrams  and  Hanratty  (1984)  to 
describe  flow  over  a  wavy  surface. 

According  to  Model  D  the  sharp  change  in  phase  angle  is  associated 
with  a  relaxation  phenomenon,  whereby  the  turbulence  does  not  respond 
immediately  to  the  oscillation  induced  variation  of  the  pressure  gradient 
Some  support  for  this  interpretation  is  obtained  from  the  measurements 
of  the  oscillation  induced  variation  of  the  mean-square  values  of  the 
turbulence  fluctuations  in  the  wall  shear  stress.  At  low  u)+  it  is  about 
in  phase  with  an  unfavorable  pressure  gradient  while  at  high  it  is 
about  in  phase  with  a  favorable  pressure  gradient  (see  Figures  5.12,  and 
5.25). 

Of  particular  interest ?  from  the  viewpoint  of  the  initial  objective 
of  these  experiments,  is  the  observation  of  a  change  of  phase  angle  over 
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a  range  of  o>+/15  approximately  equal  to  the  range  of  dimensionless  wave 
number,  a+,  for  which  similar  change  was  observed  in  the  measured  spatial 
variations  of  the  wall  shear  stress  along  a  solid  wavy  surface.  This 
suggests  that  the  same  mechanism  is  operative  and  that  the  two  experiments 
can  be  compared  by  using  the  convection  velocity  of  the  turbulence 
fluctuations  in  the  viscous  wall  region. 
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CHAPTER  7.  CONCLUSIONS  AND  RECOMMENDATIONS 

This  thesis  reports  experimental  and  numerical  results  on  how  the 
wall  shear  stress  responds  to  an  imposed  flow  oscillation  with  small 
amplitude  in  a  frequency  around  the  median  frequency  of  turbulent 
fluctuations  at  the  wall.  The  results  are  summarized  as  follows: 

a.  There  exists  a  linear  response  of  the  wall  shear  stress  to  the 
imposition  of  a  small  amplitude  flow  oscillation.  The  time 
averaged  values  of  the  wall  shear  stress  are  not  affected  by 
the  imposed  oscillation. 

b.  The  amplitude  and  the  phase  of  the  wall  shear  stress,  which 
characterizes  the  wall  response,  are  found  to  be  correlated  by 
a  single  parameter  u>+,  which  combines  the  effects  of  frequency 
and  Reynolds  number  in  the  range  of  variables  studied. 

c.  Relaxation  of  turbulence  to  the  variation  of  pressure  gradient 
in  the  wall  region  is  found  similar  to  what  has  been  observed 
for  flow  over  a  wavy  surface,  whereby  the  phase  angle  which 
characterizes  the  temporal  variation  of  the  wall  shear  stress 
undergoes  a  rapid  change  over  a  rather  narrow  range  of  a)+. 

Model  D  which  takes  account  of  this  effect  with  two  parameters 
k^  =  -25,  k^  =  3500  fits  the  data  very  well  in  this  range. 

d.  When  u)+  is  close  to  or  greater  than  the  median  frequency  of  the 
wall  turbulence,  the  measured  amplitude  of  wall  shear  stress  is 
definitely  lower  than  predicted  by  all  the  models.  There  appears 
to  be  an  interaction  between  the  imposed  flow  oscillation  and 
turbulent  fluctuations  in  the  viscous  sublayer,  which  is  not 
described  by  present  theories  of  turbulence. 
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The  present  work  is  a  preliminary  step  toward  the  understanding  of 
the  interaction  between  a  turbulent  flow  and  an  imposition  of  flow 
oscillation.  Several  interesting  problmes  emerge  from  this  work  for 
further  studies: 

a.  In  order  to  get  a  physical  understanding  of  the  relaxation 
phenomenon  it  is  necessary  to  model  the  flow  in  a  more  fundamental 
way  such  as  suggested  by  Nikolaides  and  Hanratty  (1983) . 

b.  The  mechanism  leading  to  the  lower  amplitude  of  wall  shear  stress 
at  large  to  needs  to  be  explored  both  experimentally  and 
theoretically. 

c.  At  high  and  at  larger  amplitudes  than  used  in  this  work,  the 
determination  of  the  effect  of  increasing  amplitude  on  the  flow 

is  of  great  interest,  since  there  might  exist  nonlinear  interaction 
between  turbulence  and  flow  oscillation. 

d.  The  determination  of  how  the  flow  field  can  be  affected  by  an 
imposed  disturbance  is  of  practical  importance.  Nonsinusoidal 
flow  oscillations  could  be  a  means  of  achieving  this  without 
going  to  unduly  large  flow  oscillations. 
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APPENDIX  II.  Derivation  of  Finite  Difference  Methods  for 
Variable  Grids  and  Computer  Programs. 


(A)  Finite  difference  of  variable  grids. 

A  function  of  U(y)  is  approximated  by  a  second  order  polynomial  as 


<  ,r  < 


u(y)  =  +  y  +  y  >  yj-i *  y *  yi  • 

If  the  values  of  the  function  U(y)  are  known  at  y^_p  Yj  anc*  the 

coefficients  a_.  ,  by  c_.  can  be  determined  by  solving  following  system 
equations : 


(ii. i) 


1 

yi-l 

2 

Vl 

a . 

J 

U3-l 

1 

yi 

2 

b. 

J 

= 

u. 

3 

I — 1 

yi 

yi+l 

c . 

J 

Vi 

(II. 2) 


l.  e. 


a .  = 


yJ Zi±L 


u.  ,  + 


yj-i  yi+i 


u. 


j  -  (yj.!  -  yjHyj.!  -  yJ+1>  J-i  (ya  -  -  yj+i>  J 


yi  y.i-l 


u. 


j+1 


b. 

J 


yi  +  yi+l 


u. 


yi-l  +  yi-KL 


(yj-i  -  yj)(yj-i  -  yd+i)  J-1  ^  -  y1+i>  i 


y.  +  y.  . 
i  .1-1 


U. 


(y3+l  '  yj)<y3+l  -  yJ-l}  1+1 
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c .  = 


j  <yj-i  -  V(yj-i '  yj+i>  “j'1  +  (yj  -  yj-i)(yj  -  yJ+i) 


u. 


+ 


(yj+l  ■  V(yj+1  -  *1-1>  Uj+1 


(II. 3) 


From  (3.55)  the  differences  between  two  adjacent  grids  are 

-  Vi '  4Yi  hJ'2 

Yj-1  -  yj+l  -  -  AY1  H^2  (H  +  1) 

Vl  -  ■  4Y1  hJ'1 


(II. 4) 


The  first  and  second  derivatives  of  U(Y)  are  then  approximated  as 


3U 


=  b.  +  2C.Y. 
3  ]  ] 


U 


.1+1 


_H(H  +  1) 


(H  -  1)  U.  HU. 

+ - 3 - i=L 


H 


H  +  1 


/ 


AY1*H- 


1-2 


(II. 5) 


and 


r20j±L 


_  H(H  +  1) 


H+l 


AY1  H^~2 


(II. 6) 


The  computer  programs  calculating  oscillating  flow  fields  from 


Models  A  to  D  are  as  follows: 


n  o 
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(B)  Model  A 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


10 


40 


80 

70 


this  program  solves  qu a s i -  1  am i nar  pulsating  pipe  flow  using 
cr onk-n i c  1  son  implicit  scheme, 
input:  re  reynolds  number 
vis  viscosity 

fr  frequency  of  oscillation 
ns  time  steps  per  period 
nprds  periods  to  be  proceeded 

bl  ratio  of  first  grid  size  to  the  thickness  of 
stokes  layer. 

hk  factor  of  geometric  series  of 


g  1 3  6) 
gl3. 6) 
g  1 3 .  S  ) 

—  on 
u  at  y ( j  ) 


y 

y 

y< 


output 

f or003. dat  (x, 2ei3. 5) 
f or02l . dat  —  f or028. 

velocity  profile 
f or071 . dat ( x, gi3. 6, x, 
f or072. dat ( x, g 13. 6, x, 
f or073. dat ( x, g 13. 6, x, 
for051.dat  <  x , 2g 13 . 6 ) 
time  variation  of 
program  moda 

dimension  u(200) ,  y < 200 > . a  1 f a < 200 ) 
dimension  uf < 33, 200 ) , uf f <  33,  2  > 
pi=4. #atan( 1 . ) 

print#,  'type  in  r e,  v i s ,  f r ,  ns / per i od ,  nprd s 
read#, re, vis, fr, ns, nprds 
print#,  'enter  radius  of  pipe' 
r ead# , r 
w=2 . #p i#f  r 

ymax=  099436#re## . 375 
u  s  t  a  r  =  yma  x / r 
print#  'enter  bt, hk' 
read#, bt  hk 

hl=sqrt(2. #vis/w)#ustar/bt 
m=alogiO(ymax#(hk-l)/hl+l)/aloglO(hk)+l 
print#,  '  yma x=  ' , yma x ,  'jmax=',ni 
ymax=hl*(hk#*(ni-l)-l>/(hk-l) 
n  l  =  n 1 /ymax 
initial  condition 
do  10  i=l , n i 
y  < i ) =h 1# ( nk## ( 
s  t  =  r#  s  qr  t ( w/ v  i 
u  < 1 >=0. 
do  40  i =2 , n i — 1 
tmp=-st#y ( i  ) /sqrt (2.  > 

u(  i  ) =-s  qr t ( 1  /Cl. -y( i ) ) )#exp(tmp)#sin(tmp) 
u(ni)=u(ni-l) 
print#, u  < 1 ) , u ( 2 ) , u ( 3 ) 
dt=2  #pi/ns 
l=ns/3£ 
r=ymax/ustar 
s  t  =  sqr t ( w/ v i s )#r 


grid  size 


wall  shear  stress  variation 
dat  ( x  2el3. 6) 
at  8  different  phases 


amp  1  i  tude 

phase 

zavg 


bata ( 200 ) , duy(33) 


i-U-1 )  /  (  hk-l ) 
s  ) 


u  (  4  >  ,  u  (  5  ) 


calculate  velocity  profile  at  new  time  level 
do  50  mt  =  1 . nprds 
do  60  k=l,32 
do  70  1  =  1,  1 

ip=(mt-l )#ns  +  (k-i  )#l  +  i 
a l f  a  c l ) =0 . 

Data! 1 >=0. 
call 
u  (  n  i  > 
do  80 
n=n  i  - 

u(n)=alfa(n)#u(n+l)+bata(n) 
u ( 1 )=0. 
continue 


tri  (2,  ni-l,  hi,  hk,  y,r,  st,  dt.  alfa,  bata,  ip,  u) 
=  bat a  < ni-1 )/( l. -alfa( ni-1 )) 
j=l, ni-2 
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95 

90 

60 

50 


100 
1  10 


255 

C 

c 

451 


195 


1  90 


136 


137 


510 
51  1 
■dOO 
50  1 


ne. npr d s  > 
do  95  j  =  l, 2 
u  f  f  <  k  ,  j  >  =u ( j  + 1  > 
do  90  j= 1 , n i 
u  f ( k ,  j  > =u <  j  ) 
continue 
conti nue 


go  to  60 


calculate  shear  stress  variation  at  the  wall 
do  100  k  =  l , 32 

d  u  v  <  k  >  =2 . *  u  f  f ( k , 1 ) / ( hl#st#st  > 
dea=k* i l . 25 

wrlte<3.110)  deg,duy(k) 

format(x,2el3.5) 

duv ( 33  >=duy  < l ) 

call  datf  i  t ( duy, zavg, amag, angd ) 
a  ng  d=90  +angd 
write(3,255)  y(l),amag 
write (3, 255)  y(l),angd 
wri te(3, 255)  y(i),zavg 
format! x, gl3. 6, x , g  1 3 . 6  ) 

calculate  amplitude  and  phase  profile  of  oscillating 
component  of  velocity, 
do  451  i =1 , 32 
u  (  i  ) =  uf ( i , n  i  ) 
u (33)=u( 1  ) 

call  datf i t ( u,  zavg, amag, angd  ) 

temp=angd 

do  190  1=2, n i 

do  195  1=1,32 

u ( i ) =uf ( i , j  ) 

u ( 33  >=u (  1  ) 

call  d a t f  i  t < u , z a vg , amag , ang d ) 

angd=anad-temp 

wr  i  t e ( 7 1 , 255 >  y(j),amag 

wr  l  te(72,  255)  ylVangd 

wr  i  t e < 73, 255 >  y(j),zavg 

continue 

print*,  'do  you  want  u  profile  in  a  phases  ?  0  yes,  1  no 
read*, ncr 

i  f < ncr . ne. 0 )  go  to  501 

print*,  'enter  y/r  (the  y  where  the  profile  ends)' 
r  e  a  d  *  ,  y  s 

Jr1 

iMy(j).gt.ys)  go  to  187 

j  =  j  + 1 

go  to  136 
n  k  =  j 

do  300  i=l,3 
i  i  =  i*4 
n=20*  i 

do  511  j  =  1 , n k 

write(n,510>  y ( j  ) , u f (  i  i , j  ) 

format  (  x,  g  13.  b,  x,  gl3.  6) 

continue 

continue 

continue 


print*,  'do  you  want  time  variation  of  u  at  v<i)  7  0  ves 
read* , nt  1  1 

i  f ( n  1 1  ne  0)  go  to  123 
n  =  50 

120  orint*,  'which  y/r?  type  in  j' 
read*, j 
n  =  n+  l 

do  130  k  =  1 , 32 
u(k)=uf(k,  j) 
u  <  k:  +  l  )  =  u  (  rC  1 


1  3  0 

w r  1 1 e ( n ,  1 3 1  >  dea.u(k  +  i) 

131 

format! x . 2  g 1 3 . 6) 

pr  i  nt* .  'do  on  ^  yes 

0,  no 

read* . nts  t 

if(ntst.eq.O)  go  to 

120 

123 

continue 

stop 

end 

J 


1  no ' 


on  ononnoooooooonnoooonnoon 
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(C)  Model  B 


this  program  solves  pulsating  flow  field  by  using  model  b 
crank-ni co  l  son  implicit  scheme  is  used 


input. re  reynolds  number, 

vis  vi  scos  i  ty 

fr  frequency, 

ns  time  steps  per  reriod. 

rO  radius  of  pipe, 

bt  ratio  of  first  grid  size 

stokes  layer. 


to  the  thickness  of 


output 


f  or  021 .  da  t  <  1  x,  f  13.  3,  1  x ,  4g  1 3.  6,  i3)  f  r , wpp . amag, angd, 

temp, i tr2 

of  wall  shear  s  tress . 
f  requency . 

dimensioless  frequency(wplus/15). 
amp  1 i tude. 

phase! re  lated  to  central  velocity), 
phase  of  central  ve  loc  i  ty ( re  lated  to 
pressure  gradient), 
number  of  iteration 
dat  (  x,  g  13 .6,  x,  q  13.  6 )  y(j),amag 
amp  1  i  de  prof  i  le 

for072.dat  (  x,  g  13.  6,  x  g  13.  6 )  y(j),angd 
phase  profile 


results 

f  r 

wpp 

amag 

angc 

temp 

i  tr2 
f  or07 1 . 


program  modb 

dimension  u ( 300 ) , y  <  300 ) , a  1 f  a ( 300 ) , bat  a (300 ) , du ( 33 ) , o Id  < 130 ) 
dimension  uf< ISO, 300  >,a(300),b<300), c  <  300 ) ,  strss(lSO) 
dimension  vtl( 300 ) , vt2( 300 ) , F 1 ( 300 ) , df  1 ( 300  > , dy ( 300 ) , f  ra ( 2 ) 
dimension  v(300>, f2(300>, df 2(300), f i 1(300), f i 2(300), reyd(20) 
p  i  =4.  *  a  t  an (  1.  ) 
prints  'type  in  bata' 
read*, bt 

prints,  'type  in  re, v  i  s, f r, ns, hk  ' 
read*, re  vis, fr, ns, hk 
print*,  'type  in  radius  of  pipe  ' 
read*, rO 
r  =  r  0 

w=2. *p i*f  r 

ymax  =  .  099436*re**.  875 
u  s=ymax /rO 

hl=sqrt(2. *vis/w)*us/bt 

ni=aIoglO(ymax*(hk-l )/hl+l )/aloglO(hk)+l 
ymax=hl*(hk**(ni-l)-l)/(hk-l> 
n  l  =  h  1 /ymax 

c  initial  condition 

do  10  i  =  l , n i 

10  y< i >=hl*(nk**( i— 1 )— l )/<hk— 1 > 

st=rO*sqrt(w/vis> 
u ( 1 >=0 

do  50  i =2, n i - 1 
tmp=-s t*y ( i ) / sqr t ( 2.  ) 

50  u(i)=-sqrt(l./<l.-y(i)))*exp(tmp)*sin(tmp) 

u(ni)=u(ni-l) 
do  55  i  =2 , n i 

55  dy ( 1 >=hl*nk** ( i-2> 

hk 1  =  1 . / ( hk  +  1 .  ) 
hkk 1=1  / ( hk* ( hk+l ) ) 
d  t=2  *p  i  / ns 
1  =  n  s / 32 
r=yma  x /us 
st  =  sqrt(w/vis  )*r 

call  vt(y,us,r,mini,min2,ni,vis,vtl,vt2,fl) 
do  65  J=l, mini 
65  f 1 (  i  )=vtl ( j  ) 

do  60  j=mi nl+i, n i-1 
60  f 1 ( j  >  =v  t2 ( j  ) 


70 


f  i  1  <  j  >  =f 1 <  j  ) 
continue 
f 1 ( 1 >=0. 
f i 1 ( 1 >=0. 
f 1 ( ni >  =  f l ( n i-1 ) 
f  i  1  (  n l ) =f l ( n i  -1 ) 
do  80  j=2, n i-1 

dfi(j>=(hkkl#fl<j+i>-hk#hkl#fl<j-l)+(hk-l)#fl<J)/hk) 
z  /  d  y  <  j  ) 

80  continue 

do  90  j=2,ni-l 

dfl(j>=(drl<i)-<l.+fl(i))/(l.-y(j)))/(st#5t) 

f 1 ( j  >  =  ( 1 .  +  f 1 ( j  )  )  / (  S  t#S  t  ) 

90  continue 

uf  r=u  s»vis 

wp=w#v is/(ufr#ufr) 

wpp=wp/ 15. 

u b  =  r e# v i s / ( 2  #r ) 

a  1  f  ap=wp / 15 . 

do  160  j=2, n i -1 

tmpl=fllj)/(dy(j)#dy<j>) 

tmp2=df 1( J  > / <  2. »dy  < j)  7 

a(j)=-dt#nkkl*<tmpl+tmp2) 

b(j)=l  +<tmpl-tmp2#(hk-l))#dt/hk 

c( j  )=dt#hkl#(-tmpl+hk#tmp2) 

160  continue 

i  t  r2= 1 

do  171  1=1, ns 
171  o  ld(  i  )=0. 

170  do  180  k  =  l , 32 
do  190  i  =  l ,  1 
l p=<  k-l )# l  +  i 
do  400  j=2,  ni-1 

t>ata<j)  =  (d-c(j)Hbata(j-l))/(b(j)  +  c(j)^alfa(j-l)) 

400  continue  J 

u<  ni  )=bata(ni-l )/( 1. -alfa<ni-l ) ) 
do  410  j=l, ni-2 
n  =  n  i  -  j 

410  u ( n  )  =a  I  f a ( n )#u ( n  +  l ) +ba ta < n ) 
kk  =  <  k-l )* l  +  i 
do  420  j=l , n i 
420  u  f  (  k  k , j  7=u ( j ) 

190  continue 

180  continue 

calculate  shear  stress  variation  at  the  wall  from  eamod 
d  o  430  i  =  l ,  n  s 

430  strss(i)=2.#uf<i,2)/<hl#st#st> 
do  440  i =1  ns 
def  =s  trs  s ( i ) —  o 1 d ( i ) 
if (def  gt  .OOOOtlgo  to  445 
440  continue 

do  451  i  =  1 , 32 

i i  =  i#  1 

451  du < i )=uf  < i  i  . n  i  ) 
du  < 33 )=du ( 1 ) 

call  datf i t < du . zavg, amag, angd > 
temp=angd  * 

do  450  1=1. 32 
l i  =  i#  ns /32 

450  d  u  < i )=strss( l  i) 
du ( 33 ) =du ( 1  ) 

call  datf i t ( du.zavg, amag, angd  ) 
amag=amag# 5. 02836* ( s t# rO/ r >##2/ re##  875 
anga=anga-temp 

write (21  255)  f r, wpp, amag, angd, temp, itr2 
go  to  47S 
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445  do  442  i=l , ns 
442  old( i )  =  s t rs s (  i  > 

i  t  r2= i t  r 2+1 
go  to  170 

255  f orma  t  < 1 x  f 13. 3,  1  x,  4g 13. 6, i3) 
c  profile  of  amplitude  and  phase  of  u 

47G  do  475  j=2,  n  i 
do  470  1=1,32 

i i  =  i*  1 

470  du ( i >  =u f (  i  i ,  j ) 
du ( 33  >  =du (  1  ) 

call  d a t f i t ( d u , z a v g , amag , a ng d  ) 
angd=angd-temp 
wr i te < 7 1 ,  16 1  >  y<j),amag 
wr i te < 72, 16 1 )  y(j),angd 
161  f orma t ( x , g 13. 6, x , g 13. 6 > 

475  conti nue 
stop 
end 


(D)  Model  C 


c 

r 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

c 

r 

C 

c 

c 

c 

c 

c 


this  program  solves  pulsating  flow  field 
crank-n i co  1  son  implicit  scheme  is  used. 


by  using  model  c 


input 


output 


re  reynolds  number, 

vis  viscosity, 

fr  frequency, 

ns  time  steps  per  period 

rO  radius  of  pipe, 

bt  ratio  of  first  grid  size 

of  stokes  layer. 


to  the  thi ckness 


f  or021 .  dat 


resu  1  ts 
f  r 

wpp 

amag 
anga 
temp 


(x, f 13. 3, x, 4g 1 3 . 6. i3) 
f  r,  wpp  amag, angd, temp, itr2 
wa  1  1  shear  stress. 


frequency(wplus/15). 


pressure 
itr2  num 


of 

f  requencv . 
dimesionless 
amplitude 

pha s e < r e  1  a t i ve  to  central  velocity), 
phase  of  central  ve  loci ty ( re lat i ve  to 
rad i ent ) . 
er  of  i  terat ion. 


the 


program  mode 

dimension  u( 300 ) , y ( 300 ) , a  1 f  a  <  300 ) , bat  a ( 300 ) , du ( 33 ) , o 1 d ( 130 ) 
dimension  uf( 160, 300 ) , a ( 300 ) , b ( 300 ) , c ( 300 ) , s  trss  < 160 ) 
dimension  vtl<  300 ) , v  t2<  300 ) , f 1 < 300 ) , df 1 ( 300 ) , dy  <  300 ) , f  ra ( 2 ) 
dimension  v(300>, f  2<  300  ) , df  2<  300 ) , f il( 300  > , f i 2( 300 ) , re yd ( 20  > 
pi =4. *atan<  l.  )  '  7 

pr  l  nt*  ' type  in  k  1  ' 
r ead* ,  k  1 

pr i nt^,  '  type  i  n  bata  ' 
read*, bt 

print*,  'type  in  re, v i s, f r, ns/perod, ert ' 
read*, re  vis, fr, ns, ert 

print*,  'type  in  radius  of  pipe  and  ampl.  of  velocity' 
read* , rO, gama 
r  =  rO 

w=2. *p i*f  r 

ymax  =  . 099436*re**.  875 
us=ymax/rO 

hl=sqrt(2.*vis/w)*us/bt 
hk=l  05 

ni=aloglO(ymax*(hk-l)/hl+l)/aloglO<hk)+l 
ymax=hl* ( hk** ( n i-1 )-l ) / ( hk-1 ) 
n 1  =  h 1 /ymax 
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c 

10 


50 

55 


60 


70 


80 


90 


160 

171 
1  70 


initial  condition 
do  10  i=l , n i 

y (  i  >  =  hl* ( hk## <i-l)-l)/(hk-l) 
st=rO#sqrt(w/vis) 
u  <  1  )=0 

do  50  i =2, n i -  1 
tmp=-st#y ( i ) /sqrt ( 2. ) 

u<  i  )=-sqrt  \  1.  /  (  1.  -y<  i  )  )  )#exp(.tmp)#sin(  tmp) 

u(ni)=u(ni-l) 

do  55  i =2, n i 

dy ( i ) =h 1#  nk## ( i -2 ) 

hk 1  =  1 .  / <  hk  +  1 .  > 

hkk 1=1 . / ( nk# ( hk+1 ) ) 

dt=2. #  p l / ns 

1  -  n  5  /  3  2 


r=ymax /us 
st  =  sqrt ( w/v  i  s  >#r 
call  vt <  yt  us, r, mini 


n  i 


m  i  n2, ni,vis,vtl,vt2,fl) 


x 1#  x  1 ) 


do  60  j=mi  n2+l .  .. . 
f  1  < j  )=2#vt2< ni-1 ) 
do  70  j=2, ni-1 
YP-y ( j  >*us#r 
tmp 1 =ex  p  < -y p /26 .  ) 
tmp=( 1. -tmpl ) 
x  1  =  4#y ( j )#r#tmp 
duy  =  vt 1 ( j )# v i 5 / ( : 

bd=2.#xl#(duy#duy)#tmpl#.4#y<j)#r#yp/26. 
tmp=bb/(2.#(i. -y( j ) )#<us#vis)*#2> 
f 2 ( j )=bb#2. #r#r/(vis#vis#gama#re#( 1. -tmp)  ) 
f il(J)=f  1(  j) 
f i 2 <  j ) =f 2<  j  ) 
continue 
f 1 ( 1 )=0 
f 2( 1 )=0. 
f i 1 ( 1 )=0. 
f i2( 1 )=0 
fl(ni )=fl(ni-l) 
f  2 ( ni ) =f  2  < ni-1  ) 
f i 1 ( ni )=f 1 ( ni-1 > 

Fi2(ni)=fl(ni-1) 
do  80  j=2  n i - 1 

dfl',  j)  =  <hkkl#fl<j  +  l)-hk#hkl#fl<j-i)  +  <hk-l)#fl<j)/hk) 
/dy <  J  ) 

df2(j)=(hkkl#f2<j+l)-hk#hkl#f2(J-l)+(hk-l)#f2<j)/hk) 

cont i nue 
do  90  j=2,ni-l 

f2(j)  =  (df2<j)-f2(j)/(l.  -y<j) > ) /  (  s  t#  s  t ) 
dF  l  ( j  >  =  (  d  f  l(  j  )-<  1.  +f  l(j))/(l.  -y(.j)  )  )/( 

Fl ( J )=( l 
continue 


+F1(J) )/(ST#S i) 


■y  (  j  ) )  )  /  (  s  t#  s  t ) 


uF  r  =  us*v i s 

wp=w* v is/(uFr#uFr) 

wp  p=wp / 1 5 . 

u  b  =  r e* v i s / ( 2. # r ) 


a  1 F  ap=wp/ 15. 
pc=-qama#wp#ub/uFr 
do  160  j-2.  ni-1 
tmpl=Fl(j)/(dy(j)#dy<j>) 
tmp2=dFl(j)/(2.#dy(j>) 

■'  ) =-d t# h k k 1# ( tmp 1 + tmp2 ) 

)=1. +( tmpl-tmp2#(hk-l) 


a  ( 
b  ( 
c< 

continue 


)=dt*hkl*<-tmpl+hk#tmp2> 


)#dt/hk 


i t  r 2= 1 

do  171  i  =  l , n  s 
o id ( i )=0. 
do  180  k= 1  32 
do  190  i  =  l,  1 
i p=( k-1 ># 1+i 


140 


400 


410 


420 

190 

180 


430 


440 


45  1 


450 


445 

442 

255 
47  G 

470 


16  1 
475 


do  400  j=2, n i -1 

d  =  d t*  cos ( (ip-  5)*dt)-a<j)*u<j  +  l)  +  <2.-b(J))*u(j)-c<j)*u(j-l) 

-  5*dt»f2<j>*kl*pc*cos<(ip-5)#dt> 

alfa<j)=-aij)/<b!j)+c(j)#alfa(j-l)) 

bata(j>=<d-c(j)*bata(j-l))/<b<j>+c(j)#alfa<j-l)) 

continue 

u(ni )=bata(ni-l )/( 1.  -alfa(ni-l) ) 
do  410  j  =  1 , n i — 2 
n=n i -j 

u(n)=alfa(n)*u(n+l )+bata(n) 

k  k  =  < k-1 >*  1  +  i 

do  420  i  =  l , n i 

uf  <  k k , J )=u  <  j ) 

continue 

continue 

calculate  shear  stress  variation  at  the  wall  from  eqmod 
d  o  430  l  =  1 , n  s 

strss( i )=2. *uf  < i, 2)/<hi#st#st> 

do  440  i  =  l , ns 

def=strss( i )  -  o  1  d  (  i  ) 

if(def.gt.  OOOOllgo  to  445 

continue 

do  451  i  =  l, 32 

i i  =  i*  1 

du ( l ) =u  f ( i  i , n  i  ) 
du ( 33 )=du  <  1  ) 

call  datf i t ( du, zavg, amag, angd  ) 

temp=angd 

do  450  1=1,32 

i i =i*ns/32 

du ( i )  =  s  trss (  i  i  > 

du ( 33 )=du (  1  ) 

call  datfit(duzavg, amag, angd) 
amaq=amaq*5  02836# ( s t#ro/ r  >*#2/ r e## . 875 
ango=anad-temp 

wr 1 1  e  <  2 1  255 )  fr, wpp, amag, angd, temp,  itr2 
go  to  476 

do  442  i=l, ns 
old< i >  =s  t  r  s  s (  l  ) 
l t  r2= i t  r2+ 1 
q o  to  170 

format! lx, f 13. 3, 1 x , 4  g 1 3 . 6 ,  i3) 
orofile  at  amolitude  and  phase  o-f  u 
dc  475  J  =2 ,  hi 
do  470  1=1, 32 
l  i  =  i#  1 

du(i)=uf(ii,j) 
du ( 33 )=du  <  1  ) 

call  datf i t ( du, zavg, amag, angd  ) 

anad=angd-temp 

write! 71, 161 >  y(j),amag 

wr i t e ( 72, 16 1  >  y!j),anqo 

f  orma t ! x , gl3. 6, x, gl3. 6) 

continue 

stop 

end 


no  o  o 
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(E)  Model  D 


c 


c 

r 

C 

w 

c 

c 

c 

c 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


10 


50 

55 


this  program 
P  y  model  c( 
input  re 
v  1  s 
f  r 
n  s 

nprds 

bt 


solves  turbulent  pulsating  flow  in  a  pipe 


hk 
k  1 
k  1 

output. f  or 002 
f  or003 
f  or004. 
f  or021 
f  or022 
f or023. 


to  the  thickness 


i  (  J  ) 


in  time 


reynolds  number 
viscosity 
f  requency 

time  steps  per  period 
periods  to  be  proceeded 
ratio  of  first  grid  size 
of  stokes  layer. 

factor  of  geometric  series  of  grid  size 
constant  in  model  d 
relaxation  constant  in  model  d 
dat<x,  i  3 , 4g  1 3 . S  )  j,y<j>,vtl(j),vt2(j).f 
e  d  d  y  v  i  s  c  o  s  i  t  i  e  s 
dat < x , 2g 13. S >  deg, stress 

shear  stress  variaation  at  the  wall 
dat ( x , 2q  13. 6  )  deg,p(eff) 

effective  pressure  gradient  variation 
dat(x4gl3  6)  y(j),mean,amag,angd 
profile  of  mean  ,ampl  and  phase 
dar,  (x2gl3.6)  y(j),u(ii,j) 
profile  of  u  in  s  phases 
dat < x , 2g 13. 6  )  deg  u(ii,j) 
u  variation  in  time  at  y(j) 
f o r024 . d a t < x  2g 13 . 6 )  y(j),b(j)  (re.  stress) 
profile  of  reynolds  stress 

program  mod2 

dimension  u<  300  ) ,  y  < 300 ) , a  1 f  a ( 300 ) . ba ta ( 300  > , du (33 ) 
dimension  uf (  ISO, 300 ) , a ( 300 ) . b ( 300 ) , c(300>, strss< ISO) 
dimension  vtl<  300 ) , vt2<  300 ) , f  1 ( 300  > , df  1 ( 306 ) , dy <  300 ) 
dimension  v  (  300  > ,  f  2( 300  ) ,  df 2( 300  )  ,«*  i  1  <  300  )  ,  f  i2(  300  ) 
pi =4. #a  tan (  l.  ) 

re, vi  s, fr, ns/period, nprds  ' 
ns , nprds 

radius  of  the  pipe  and  amplitude  of  velocity 


pr  i  n  t# ,  '  ty pe  in 
read#, re  vis, fr 
print#,  '  type  i  n 
r ead# , rO  g  ama 
print#,  'type  in 
read# , cr t 
r-r  0 

w  =  2  #  p  l  #  f  r 

ymax  =  099436#re## 

us=ymax/r 


converge  cr  i  ter i a 


975 


select  grid  size 
pr  i  nt#,  ’bata, hk  ' 
read#, bt, hk 
aa=26. 
c  k  —  A 

hi=sqrt(2. #vis/w)#us/bt 


ni=aioglO(ymax#(hk-l)/hl+l)/aloglO(hk)+l 
print#,  'ymax=  ' , ymax,  '  jmax=  ' , n i , ’ h 1=  ' , hi 
vmax=hl# ( hk## ( n i-l )-l ) / ( hk-1 ) 
h  1  =  h 1 /yma  x 

print#,  'hl=',hl,  'ymax=',ymax 


initial  condition 
do  10  i  =  l , n i 

y( i )  =  hl#(hk##(  i  —  1  )-l >/<hk  —  1 > 
s  t  =  r#  sqr t (w/vis) 
u  (  1  )  —  0 

do  50  i =2, n i - 1 
tmp=-st#y< i ) / s qr t ( 2 .  ) 

u< i ) =-sqr t ( 1. /( 1. -y( i ) ) >#exp(tmp)#sin(tmp) 

u<  ni  )=u( ni-i )  v 

do  55  i  =2, n i 

dy (  i  )  =h l#hk## ( i-2 ) 

h k 1  =  1  / ( hk  +  l .  ) 

hkk 1  =  1 . / ( hk# ( hk  +  l )  ) 

dt=2. #pi/ns 

l=ns/32 

r  =  ymax /u  s 

s  t  =  sqr t (w/vis )#r 
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c 

c 


SO 


70 


z 


z 

80 


90 


c 

c 


105 

100 


120 

130 

140 

145 

146 

150 

151 


vt 1 . vt2, f 1 , aa . ck ) 


aa 


calculate  u  at  new  time  level 
call  v t ( y , u s , r , m l n 1 , m l n2, n i , v i s , 
do  60  j =m l n2+ 1 , n i 
f  1  <  i  >=2#vt2( ni-1 ) 
do  70  j=2, n i-1 
yp=y< j  >  *  u  s  *  r 
tmpl=exp< -yp/aa ) 
tmp= ( 1 . -tmp l > 
x  l  =  ck#y  <  j )*r#  tmp 
d  u  y =v  1 1  ( j  )*  v  i  s/ ( x  1*  x  1  > 
bb=2. *x  l*<duy*duy)*tmplwck#y<j )#r# y p/ 
tmp=bb/(2.H<l.-y<j))#(us#vis)##2> 
f  2  ( j )  =  bb#2. *r#r/<vis#vis#gama*re#< 1. -tmp) > 
f ill j)=f 1( j) 
f  i2( j >  =  f  2 ( j  ) 
continue 
f 1 ( 1 )=0 
f 2< 1 )=0 
f i 1 ( 1 )=0. 
f 12 ( 1 >=0. 
f 1 ( ni >  =  f 1 <  ni-1 ) 
f 2( ni >=f 2( ni-i  > 
f 1 1 (nl )sf 1 <  nl-1 ) 
f i 2 (  ni )=f 2<  ni-1 )  . 
do  80  j=2, ni-1 


>  =  (h(ckl#fl<j  +  i>-hk#hkl#fl(j-l)  +  <hk-l>#fl<j)/hk) 
>=<hkkl#f2<j+l)-hk#hkl#f2(j-l)+<hk-l)*f2(j)/hk) 


df  l  ( 

/dy  ( 
df  2< 

/dy  ( 

continue 
do  90  ,i=2 ,  ni-1 


do  90  j =2 , ni-1 

f 2( j  >  =  (df 2< j )-f 2( j )/( 1 . -y  Cj  >  >  > 
df  1  7  j  )  =  <  df  1  7  j  >-<  1 .  -tf  1  (  j  ) } /( 1  - 
+f  l  (  j  >  >  /  < ft*st ) 


j  >  =  <c 
f i<  j )  =  (1. 
continue 
ufrO=  099436*re** 
w 
u 


)/<  st* st ) 
y(j))>/(st#st) 


=w*  vis/(ufrO#ufrO) 
=rs«v i 5/ ( 2.  nr ) 


875*v is /r 


a  1  f  ap=wp/ 15. 


calculate  effective  pressure  gradient 

prints  ' type  in  k  1 ,  k  1  ' 

read#, k  1 , k  1 

pc=-gama*wp*ub/ufrO 

a  1 f  ap=2*a l F  ap#k  1 

pk=(klHdt)/<alfap+dt> 

vk  =  <  a  1 f  ap-d t ) / ( a I f  ap+d t  > 

v( 1 >=0. 

i  t  r=  1 

continue 

do  100  i=2, ns+1 

i  i  =  <  i  ti — 1  )*ns+i-l 

v< i  >  =  pk#  pc#  < cos( 1 i#d  t )  +  cos( < i i-1 )#d  t ) )  +  vk# v ( i-1 ) 

v ( 1 )=v  <  ns  + 1 ) 

do  120  i=l . ns 

d 1 1 a  =  ab s ( v < i ) -d f  2 ( i  )  > 

if < d 1 ta. g t . crt >  go  to  130 

continue 

go  to  145 

do  140  i=l , ns 

DF2 ( I >=V<  I  > 

i tr=l+i tr 

go  to  105 

do  146  i = 1 . 32 

i i  =  i#  1 
deg= i#  1 1 . 25 
du ( i ) =v (  i  i  ) 

write(4,l50)  deg,du(i) 
format ( x . 2gl3. 6) 
du ( 33 )=du ( I  ) 

call  datf i t < du, zavg, amag, angd  ) 
write<4,151)  zavg , amag. angd 

format!/,  x,  'mean=',gl3.6,x,  'ampl  =  ',gl3.6,  x,  'phas=',gl3.6> 
print#,  '  i tr=  ' ,  i  tr 
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160 


z 


400 


410 


420 

190 

180 

170 


430 


450 


460 


c 

r 


451 


470 


161 

475 

c 


186 


187 


do  160' j=2, n i-1 

tmpl=flij)/(dy(j)#dy(j)> 

tmp2=df 1  T  j )  /  ( 2. #dy  ( j )> 

a  <  j  ) =-d  t#  nk  k 1#  ( tmpl  +  tmp2) 

b( j  >  =  1 . +( tmpl-tmp2#(hk-l ) )  #  d  t  /  h  k 

c(  j  >=dt#hkl#(-tmpl+hk#tmp2) 

continue 

do  170  m=l , nprds 
do  180  k=l , 32 
do  190  i  =  1 ,  1 
i p=< k-i )* 1+ i 
do  400  j=2, n i-1 

d=dt#cos<<ip-.5>#dt)-a<j)#u(j+l)+(2.-b(j>)#u<j)-c<j)#u(j-l) 

-  5#dt#f2(j)#(v( i p ) +v ( ip+1) ) 

a  1  f  a( j  >=-a ( j )/<b(j>+c<j)#alfa<j-l>> 

bata<j>=<d-c<j)#bata(j-l))/<b(j>+c<j)#alfa(j-l>> 

continue 

u(ni >=bata(ni-l)/( 1.  -alfa(ni-l) ) 
do  410  j=l ,  ni-2 


n=n l - j 

u(n)=alfa(n)#u(n+l)+bata(n) 

l f ( m.  ne. nprds  )  go  to  190 

kk=<  k-1  )#  1  +  i 

do  420  J  =  1 , n  i 

uf(kk, j 1 =u ( j  ) 

conti nue 

continue 

continue 

calculate  shear  stress  variation  at  the  wall  from  mod2 
do  430  1=1, NS 

strssl i >=2. #uf  < i,2)/<hl#st#st> 
do  450  1=1,32 
I  1  =  I#L 

du ( i )  =  s t r ss (  i  i  ) 
d  u ( 33 ) =d  u (  1  ) 
do  460  i = 1 ,32 
deg= i# 1 1 . 25 

write(3,150)  deg,du(i) 

call  da t f i t ( du , za vg , amag , angd  ) 

angd=90. +angd 

w r i t e ( 3 ,  151)  zavq, amag, angd 
amag=amaq*5. 0233b# (st#r0/r)##2/re##. 875 
print#,  '  t ( 0 ) =  ' , amag 


profile  of  amplitude  and  phase  of  u 
do  451  i=l,32 
l i  =  i#  1 

du ( i )=uf ( i i , n  i  ) 
d  u ( 33  >  =d  u  <  1  > 

call  d a t f i t ( du , za vg , amag , angd  ) 

temo=angd 

do  475  j  =2  n  i 

do  470  i=l , 32 

i  i  =  i#  1 

du ( i ) =u f  <  i  i , j  ) 

continue 

du ( 33  )=du  <  1  ) 

call  da t f 1 t < du , zavg , amag , angd  ) 
arngd  =  angd-temp 

wr i te < 21 ,  161)  y(j), zavg, amag, angd 

f  orma  t  <  x , 4g 1 3 . 6  > 

continue 


print#, 'do  you  want  u  profile  in  8  phases  ?  0  yes,  l,no' 
read# , ncr 

i f ( ncr  ne. 0  )  go  to  501 

print#,  'enter  y/r  (the  y  where  the  profile  ends)' 
r ead# , y s 

j  =  1 

i f ( y ( j ) . g t . y s  )  go  to  187 
j  =  j  +  l 
go  to  186 
nk=j-l 


144 


511 

500 

501 


121 


131 

123 


530 


525 

520 

521 


do  500  1=1,8 
i i  =  i#  n  s  /S 
do  511  j=l , nk 

wr  i te<22,  150)  y ( j ) , uf (  i  1 , j  ) 

continue 

continue 

continue 


print#,  'do  you  want  time  variation  of  u  at  y<J)  ?  0  yes 
read# , n  1 1 

i f  <  nt  1 .  ne  0 )  go  to  123 
print#,  'which  y/r  ?  type  in  j' 
read#, j 
do  131  1=1,32 
i  i  =  i  #  1 


1 


deg=i# 1 1 . 25 

wr i te ( 23, 150 )  deg, uf ( i i, J ) 
print#,  'go  on  ?  0  yes  ,  l  no' 
read#, ntst 

i  f 1 ntst  eg. 0)  go  to  121 
continue 


print#,  'do  you  want  reynolds  stress  profile  at  8  phases  ?' 
pr i nt#,  '0  yes,  1  no  ' 
read#, nt2 

i f ( ntz. ne. 0 )  go  to  521 
do  520  i=l, 8 
i  i  =  i#ns/8 
do  530  j=2, nk 

a ( J  >  =  <  u  f (  ii, j  +  l)#hkki-hk#hkl#uf( ii,  j-l)  +  (hk-l)#uf(  ii,  J)/hk) 

/dy ( j  ) 

b(^)=-(fil(j)#a(j)-fi2(j)#v(ii))#gama/re 

do  525' j=l , nk 

wr  i  te < 24,  150 )  y(J),b(j) 

continue 

continue 

stop 

end 


no ' 


(F)  Subroutines 


c 

subroutine  tri(m, n,hl,hk,y,r,st,dt,alfa,bata,  i,u) 
dimension  u ( 200 ) , a l f  a  <  200) , ba  t a  <  200  ) , y  <  200 ) 
n  k  l  =  l  /  (  n  k  + 1 ) 

bkk  1  =  1 .  / ( hk# ( hk  +  l ) ) 

do  300  j=m, n 

dv=ni#hk## ( j-2 ) 

stdy=dt/ ( st#st#dy#dy ) 

stdyym=dt/(2. #st#st#dy#< 1  — y  <  j  > ) ) 

a=(-stdy+stdyym)#hkkl 

b  =  l .  +(stdy+stdyym#(hk-l))/hk 

c=(-stdy-stdyym#hk)#hkl 

d=dt#cos( ( i-1+0. 5)#dt)-a#u(j+l)+<2. -b)#u<j)-c#u(J-l) 

alfa< j  )=-a/ ( b+c#a If a( j-1) > 

bata( j  >  =  (d-c#bata(j-l) )/(b  +  c#alfa(j-l) ) 

300  conti  nue 
r  e  t  u  r  n 
end 


c 


132 


133 

134 


C 

C 


210 

I 

250 


/ 

6 

5 

c 


subroutine  vt(ytus,r,minltmin2, ni.vis.vtl.vt2.fi) 

dimension  vtl (260), vt2<206>, f i (206), y(200^ 

d ef  m 1  =  1 0000 

def m2=10000. 

do  132  j  =  2, n i -1 

vp=y ( j )»u5#r 

tmpl=exp( -yp/2S.  ) 

vt 1 ( j )  =  . 5* ( -1 . +sqrt( 1. +. 64*yp**2* (1. -tmpl )*#2#  < 1. -y(j) >  >  > 
tmp2=<2. -y< j > >#< l  +2 . # ( 1 . -y < j ) ) ##2 ) 
vt2< j  )  =  .  4*yp* tmp2/6. 
x  1  —  4#r#y  <  i  )*  (  1 .  -  tmpl) 

x  1  1  =  .  2*  r*y ( j >*tmpl*yp/<26. *  ( 1. -y<  j ) )*  <  us* v i s )**2) 

dmvy=vis*vtl(j)/(xl*xl) 

tmp  =  2  * x  1* x  1  l*dmvy**2 

l  f ( j .  g  t . 20  )  go  to  132 

f i<  j)=<2,*vtl(j)+tmp)/(l. -tmp) 

do  134  j  =20 ,  n i - 1 

def l=abs ( vt 1 <  j )-vt2<  j ) ) 

def2=abs<  f i ( j>-2*vt2ini-l ) ) 

i  f ( def 1 . qe. def ml )  go  to  133 

defml=defl 

m i n 1= j 

i  f ( def 2.  ge. def m2)  go  to  134 

def m2=def 2 

m i n2= j 

continue 

return 

end 


suorout  ine  datf  it<  z.zavgamag, angd) 
dimension  a  (  2 ,  2 ) , b  (  2  > , z ( 33 ) , x  ( 33 ) 
n  p  =  33 

pi =4  *a t an (  1 .  ) 
n  =  2 

do  210  1=1, 33 
x (  i  )  =  ix . 03 125*2* p i 
do  1  i  =  1 , n 
do  2  j= 1 , n 

a (  l  ,  j  )  =0 
b  i  l ) =0 

sumava=0. 
d o  250  i  =  1 , 32 
sumavq=sumavg+z< i ) 
z  a  v  q  =  s  u  m  a*  v  g  /  v  np-1  ) 
do  5  m=l , np 

l  l  =0 


do  S  j  1  =  1,  n,  2 
i i = i l +1 


b(ij)=b<jj)+(z(m)-zavg)*sin(x(m)*ii) 
b(  j )  + 1 )=d( jj  +  l )  +  <  z(m)-zavg)*cos( x(m)*i i ) 
mm=o 


do  7  kk  =  l,n,2 
mm=mm+ 1 

a  <  j  j ,  k k  >=a ( j  j ,  kk)+sin(x(m)*ii )*sin(x(m)*mm) 

a  (  j  j  ,  kk  + 1 >=a <  j  j ,  kk  +  1 )+sin(x(m)*ii)*cos(x(m)*mm) 

a(ji+l,kk)=a(jj+l,kk>+cos<x(m)*ii)*sin(x<m)*mm) 

a(jj+l,kk+i)=a(jj+l,kk+l)+cos(x(m)*ii)*cos(x(m)*mm) 

continue 

continue 

continue 


de 1 t=a < 1 ,  l )*a<2, 2>-a< 1, 2)*a(2, 1 ) 
deltl=b( 1 )*a  <  2, 2 )  -  b  <  2  )#  a  ( 1, 2) 
de 1 t2=b ( 2 )*a ( 1 . l)-b(l)*a(2, 1) 
b ( 1 ) =de 1 1 1  /de  1 1 
b ( 2 )  =  d  e  1  t2/de  1 1 
amag=sqrt( b( 1 >**2+b(2)**2) 
bb  =  -b ( 1  ) 
b  b  b  =  b ( 2 ) 

ang=atan2( bb, bbb ) 
angd=anq*360  / (2. *p  i  ) 
r e  tu  r  n 
end 
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APPENDIX  III.  Sampling  Programs  in  LSI-11 


(A)  Sampling  program  for  steady  flow. 


************************************************** 


PROGRAM  MCRTS 

THIS  IS  THE  FORTRAN  CONTROLLER  PROGRAM 
FOR  ANALOG  TO  DIGITAL  CONVERSION  OF 
MULTIPLE  CHANNELS  AT  SPEEDS  UP  TO 
20000  CHANNELS  PER  SECOND  THIS 
PROGRAM  SETS  THE  OUTPUT  FILE, THE  ERROR 
FLAGS, THE  DATA  ARRAY  AND  CALLS  THE 
MACRO  SUBROUTINE  MC-RTSB  .  MAC  WHICH 
STARTS  THE  CLOCK,  AND  COLLECTS  THE  DATA 


THE  OUTPUT  FILE  IS  SPECIFIED  USING  THE 
GETSTR  COMMAND,  THIS  MUST  BE  COMPILED 
USING  THE  NOSWAP  OPTION. 


THE  OUTPUT  FILE  IS  FORMATTED  15,15 
WITH  THE  FIRST  A  DATA  POINT  NUMBER 
AND  THE  SECOND,  THE  VALUE  OF  THE  DATA 
POINT  IN  2 ’ S  COMPLEMENT. 

A  SEPARATE  PROCRAM  IS  NECESSARY  TO 
CONVERT  THESE  NUMBERS  TO  VOLTAGES . 
THIS  SAVES  MEMORY  SPACE  IN  THE 
COLLECTION  PROGRAM  FOR  THE  LARGEST 
NUMBER  OF  DATA  POINTS 


THE  NUMBER  OF  CHANNELS  AND  THE  NUMBER 
OF  DATA  POINTS  ARE  READ  FROM  THE 
TERMINAL  AND  ARE  TRANSFERED  TO  THE 
SUBROUTINE. 


THE  CLOCK  SPEED, THE  DATA  AQUISITION  RATE, 
IS  SET  BY  SELECTING  THE  MULTIPLE 
OF  10000/ SEC  DESIRED. 

I .2.  RATE-  MULTIPLE /  1  0 00 0 
MULTIPLE  MUST  BE  A  POSITIVE  INTEGER 
AND  IF  FINER  INTERVAL  SPLITTING  IS  NEEDED 
THE  MACRO  PROGRAM  CAN  BE  EDITED  TO 
ALLOW  DIVSIONS  OF  10**5  OR  10**6/SEC. 


****************************************************** 


c 


c 


c 


1 

10 

c 


c 


PROGRAM  MCRTS 

INTEGER* 2  DATA (  1  2  000  ) 

DATA  DATA/ 1 2000*0/ 

LOGICALM  DATFILU2)  ,  ERR 

COMMON  /AD/  DONE , ERCLK , ERADC  ! DONE  FLAG, ERROR  FLAGS 

DONE-Q 

ERCLX=0 

ERADC=0 

WRITE (7, 10) 

FORMAT! IX , 'TYPE  IN  THE  NAME  OF  THE  OUTPUT  FILE') 

CALL  GETSTR! 5 , DATF I L , 11 , ERR ) 

IF (ERR)  GO  TO  1 


C 

C 


OPEN! UNI T^ 3 , NAME = DATF I L , TYPE= ‘ NEW ' ) 


.'WRITING  WITH 

!  DEVI CE  3  WILL 

! NOW  GO  INTO  OUTPUT  FILE 


o  o  o 


C 

C 

20 


C 

30 

C 

C 


40 

C 

C 

99 

C 


50 

C 

89 

60 

C 

79 

70 

C 


80 

00 


*  *  *  * 
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VR  I  TE  (  7  ,  2  0  ) 

FORMAT ( 1 X , ‘ TYPE  IN  THE  NUMBER  OF  CHANNELS 
AND  THE  NUMBER  OF  POINTS  PER  CHANNEL') 
READ  <  5  #  * )  NCHS , NPTS 


Z 


VR ITE ( 7 , 3  0 ) 

FORMAT ( 1 X , ' TYPE  IN  THE  MULTIPLE  OF  10000/ 
FOR  THE  SAMPLING  RATE • 


READ ( 5 , * )  MULT 


VR ITE  <  7 , 4  Q ) 


FORMAT( IX; 'DATA  AQUISITION  IN  PROGRESS') 

CALL  MCRT5B ( NPTS , NCHS , MULT , DATA ) 

IF  (DONE  EQ . 0 )  GO  TO  99 

IF  (ERCLK.LT. 1)  GO  TO  89 
VR ITE ( 7 ; 5  0 ) 

FORMAT( IX; 'CLOCK  OVERFLOW  ERROR') 

IF  ( ERADC . LT. 1 )  GO  TO  79 
VR I TE ( 7 ; 60) 

FORMAT ( IX ,  1  A/D  CONVERTER  ERROR') 

WRITE (7, 70) 

FORMAT( IX; 'DATA  COLLECTION  COMPLETE') 

DO  100  1=1 , NPTS*NCHS 

VR ITE ( 3 ; *  >  DATA(I ) 

WRITE (3/80)1, DATA  < I ) 

FORMAT ( 15 ; I  5 ) 

CONTINUE 


STOP 

END 


SUBROUTINE  MCRTSB . MAC 

THIS  SUBROUTINE  IS  CALLED  TROM  MCRTS  FOR 
CALL  MCRTSB (NPTS, NCHS, MULT, DATA) 


WHERE: 


NPTS  =  NUMBER  OF  DATA  POINTS 
NCHS*  NUMBER  OF  CHANNELS 
MULT=  MULTIPLE  OF  10000/SEC 
DATA=  DATA  ARRAY 


PER  CHANNEL 


SiSlflSS  4  15  USED  T0  POINT  TO  THE  DATA  ARRAY 
REGISTER  5  IS  THE  POINTER  FOR  THE  ARGUMENTS  IN 
THE  FORTRAN  CALL 

THE  OTHER  REGISTERS  ARE  NOT  USED 


**************************„** 


*********************,  *****»., 


. SBTTL  MCRTSB. MAC 
. GLOBL  MCRTSB 
.GLOBL  COLOOP 


AD3TAT= 177000 
ADBUF  =ADSTAT+2 


I  THE  ADDRESS  OT  THE  A/D  BOARD 
;THE  ADDRESS  OF  THE  DATA  REGISTER 


CLSTAT= 170420 
CLBUF  =  C  LSTAT  +  2 


) THE  ADDRESS  OF  THE  CLOCK 
) THE  ADDRESS  OF  THE  CLOCK  COUNTER 


C  LVECis  4  4  0 
PR7  =  3  4  0 


.THE  CLOCK  INTERUPT  VECTOR 
; INTERUPT  PRIORITY  IS  7 
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DONE 
ERCLK  * 
ERADC 


DATA 

NCHS 

NPTS 

MULT 


MCRTSB : 


COLOOP: 


1$  : 


22$  : 
2$  : 


3$  : 


EXIT: 


10$  : 


DEVPR I r 5  ; DEV I C E  PRIORITY  IS  5 

PSECT  AD , RW , D , GBL , REL , OVR 
WORD  0 
WORD  0 

•WORD  0  ; SET  UP  FLAGS  IN  /AD/  COMMON 

.PSECT  $DATA  RV , D , LC L , REL , CON 
.WORD  0 
.WORD  0 
WORD  0 
.WORD  0 


.PSECT  $CODE  RW, I , LCL , REL , CON 

MCALL  . INTEN  ; LOWERS  PRIORITY  ON  INTERUPT 


MOV  9<R5)+,R0 
MOV  9  <  R5 ) ♦  , NPTS 
MOV  0(R5)+,NCH5 
MOV  9 ( R5 ) f i MULT 
MOV  (R5  )  + ,DATA 


IGNORE  NUMBER  OF  ARGUMENTS 
GET  NUMBER  OF  POINTS 
GET  NUMBER  OF  CHANNELS 
GET  MULTIPLE  FOR  CLOCK 
GET  ADDRESS  OF  DATA  ARRAY 


MOV  ICOLOOP ,  0#CLVEC  ;  LOAD  INTERUPT  LOCATION 

MOV  #PR7 , 0#CLVEC  +  2  ; AND  PRIORITY 

CLR  @#ADSTAT 
TST  9# ADBUF 


CLR  0#CLSTAT 
NEG  MULT 
MOV  MULT , 9#CLBUF 
MOV  #133 ,0#CLSTAT 


INITIALIZE  CLOCK 
CHANGE  MULT  TO  2‘S  COMPLEMENT 
MOV  MULT  INTO  CLOCK  COUNTER 
START  CLOCK  10000/SEC 
INTERUPT  ENABLED 


RTS  PC 


; RETURN  TO  MAIN  PROGRAM  TO 


, WA IT  FOR  INTERUPT 


INTEN  DEVPR I 
3IC  #200 ,0#CLSTAT 
MOV  NCHS , R5 


LOVER  DEVICE  PRIORITY 
CLEAR  CLOCK  OVERFLOW  BIT 
SAVE  NUMBER  OF  CHANNELS 


MOV  DATA , R4 
MOV  # 1 , 9# ADSTAT 
TSTB  9#ADSTAT 
BPL  1$ 

MOV  §#ADBUF , (R4)+ 
DEC  NCHS 
BEQ  3  % 

ADD  #401, 9#ADSTAT 
TSTB  9#ADSTAT 
BPL  2$ 

MOV  9#ADBUF , <R4)  + 
DEC  NCHS 
BNE  22$ 

DEC  NPTS 
BEQ  EXIT 
MOV  R4 , DATA 
MOV  R5 , NCHS 
RTS  PC 


; MO VE  DATA  ADDRESS  TO  R4 
, LOAD  CHANNEL  0  AND  START  CONVERSION 
; VA IT  UNTIL  COMPLETE 

; SAVE  DATA  POINT 


i SET  SEQUENTIAL  MODE 
; TEST  IF  DONE 

; SAVE  NEXT  DATA  POINT 
; DECREMENT  NCHS 
; ANY  MORE  CHANNELS 
; NO  DECREMENT  POINTS 
; EX  IT  IF  DONE 

SAVT  DATA 

•SAVE  CHANNEL  COUNT 
; RETURN  FOR  MORE  DATA 


CLRB  9#CLSTAT 
BIT  #10000 ,9#CLSTAT 
BEQ  10$ 

MOV  $2 , ERCLK 
TST  9#ADSTAT 
BPL  RETURN 
MOV  #2, ERADC 
MOV  #1 , DONE 
RTS  PC 
.  END 


STOP  CLOCK 

CHECK  CLOCK  ERROR  BIT 
NO  SET  CONTINUE 
SET  CLOCK  ERROR 
AD  ERROR  BIT  SET  ? 

NO  RETURN 
SET  A/D  ERROR 
SET  DONE 


RETURN: 


UuUuUUOu 


(B)  Sampling  program  for  time  interval  between  two  trigger  signals 


c************************************************************* 

C  PROGRAM  TTMF . FOR 

THIS  IS  A  FORTRAN  DRIVER  FOR  MEASURING  MEAN  TIME 
INTERVAL  BETWEEN  TVQ  EXTERNAL  TRIGGER  SIGNALS, 

NPTS  NUMBER  OF  PERIODS  TO  BE  MEASURED 

NN  NUMBER  OF  POINTS  TO  BE  USED  IN  DATCL , FOR 

FOR  A  PERIOD. 

MULT  SAMPLINIG  RATE  TO  BE  USED  IN  DATCL. FOR. 

BEFORE  RUN  THIS  PROGRAM , LINK  IT  WITH  TTM.MAC 
************************************************************** 

PROGRAM  TTMF 
INTEGER*  2  DATA(5Q0 ) 

DATA  DATA / 500*0/ 

LOGICAL* 1  DATF I L ( 1 2 ) , ERR 
COMMON  /  AD / DONE , ERCLK 

DONE=0 
ERC  LX= 0 
C 

1  WRITE (7,10) 

10  FORMAT( IX , ‘TYPE  IN  THE  NAME  OF  THE  OUTPUT  FILE1) 

C 

CALL  GETSTR  <  5 , DATF I L , 1 1 , ERR ) 

IF (ERR)  GO  TO  1 

OPEN(  UN  IT-  3  ,  NAME = DATF  I L , TYFE= 1  NEW  *  ) 

C 

WRITE (7, 20) 

20  FORMAT< IX , 'TYPE  THE  NUMBER  QF  PERIODS  AND  POINTS  PER  PERIOD’ > 
READ ( 5 , * )  NPTS ,NN 
WRITE (7 , 30 ) 

30  FORMATC IX, ’PERIOD  MEASUREMENT  IN  PROGRESS') 

CALL  TTM( DATA, NPTS) 

9?  IF (DONE . EQ. 0 )  GO  TO  99 
I  F ( ERCKL . LT . 1)  GO  TO  89 
WRITE (7 , 40) 

40  FGRMAT( IX, 'CLOCK  OVERFLOW  ERROR') 

39  WR ITE ( 7 , 5  0 ) 

50  FORMAT( IX , 'PERIOD  MEASUREMERNT  IS  OVER’) 

SUM-  0 

DO  100  1=1, NPTS 
100  SUM= SUM + DATA ( I  ) 

COUNT= SUM /NPTS 
MU LT= COUNT /NN 
STDV=0 

DO  200  1=1, NPTS 

200  5TDV=STDV+ < DATA ( I ) -COUNT) **2 
STDV = SORT ( STD V /NPTS ) /COUNT 
WRITE (3,*)  COUNT, STDV, MULT 
C  110  FORMAT ( 1 X ,  * COUNT = ',16,' STDV= ’  , F 3  .  4  ) 

STOP 

END 


ft********************************************************************* 
SUBROUTINE  TTM.MAC 

THIS  SUBROUTINE  IS  CALLED  FROM  TTMF. FOR. 

CALL  TTM( DATA, NPTS) 

WHERE  DATA=DATA  ARRAY. 

NPTS  =  P ER IODS  TO  BE  AVERAGED 

***************************************************t***************** 

. SBTTL  TTM.MAC 
.  GLOBL  TTM 


CLSTAT= 1 70420 
CLBUF  =  CLSTAT+2 


. THE  ADDRESS  OF  THE  CLOCK 

;THE  ADDRESS  OF  THE  CLOCK  COUNTER 


CL VEC  =  4  4  4 
PR7  =  3  4  0 


iTHE  CLOCK  INTERRUPT  VECTOR 
; INTERRUPT  PRIORITY  IS  7 
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DONE: 
ERCLX : 


DATA 

NPTS 

FLAG 


TTM* 


CLKISR: 


1  $  * 


EXIT 


RETURN: 


DEVPR I  =  5  ; DEVICE  PRIORITY  IS  5 

.PSECT  AD , RV , D , GB L / RE L , 0 VR 
.WORD  0 
WORD  0 

.PSECT  SDATA  RV # D # LCL , REL , CON 
WORD  0 
.WORD  0 
.WORD  0 

.PSECT  $CODE  RV , I , LCL , REL , CON 
. MCALL  . INTEN 


MOV  § ( R5 ) + , RO 
MOV  <  R5 )  +  , DATA 
MOV  0 ( R 5 ) +  r NPTS 
MOV  #0 , FLAG 
MOV  #CLKISR,@#CLVEC 
MOV  IPR7  , (?#CLVEC  +  2 


IGNORE  NUMBER  OF  ARGUMENT 
GET  ADDRESS  OF  DATA 
GET  NUMBER  OF  POINTS 
SET  FLAG  ZERO 
LOAD  INTERRUPT  LOCATION 
AND  PRIORITY 


CLR  9#CLSTAT 

MOV  #60036 ,9#CLSTAT 

RTS  PC 


; INITIALIZE  CLOCK 
; SET  CLOCK, ST2  GO  ENABLE , 
i INTERRUPT  ON  ST2 , 1 OKHZ , MODE  3 
; RETURN  TO  MAIN  PROGRAM  TO  VAIT 
; FOR  INTERRUPT 


. INTEN  DEVPR I 
TST  FLAG 
BNE  1? 

MOV  #1 , FLAG 

BIC  #  1  00000  ,@#CLSTAT 

RTS  PC 


; LOVER  DEVICE  PRIORITY 
t TE ST  FLAG 

; IGNORE  FIRST  5T2  EVENT 
; SET  FLAG  ONE 
;  CLEAR  ST2  FLAG 

; RETURN  TO  MAIN  PROGRAM  TO  VAIT 
; FOR  NEXT  ST2  EVENT 


MOV  DATA , RO 
MOV  SiCLBUF, ( RG ) * 

MOV  RO , DATA 

BIC  #100000 ,8#CLSTAT 

DEC  NPTS 

BEQ  EXIT 

RTS  PC 


, SEND  COUNTS  TO  DATA  ARRAY 

;  CLEAR  ST2  FLAG 

; DECREMENT  NUMBER  OF  POINTS 

; NO  MORE  POINT,GO  TO  EXIT 

iMORE  POINTS , RETURN  TO  MAIN  PROGRAM 

; TO  VAIT  FOR  NEXT  ST 2  EVENT 


BIC  #20000, 0#CLSTAT 
BIC  1 1  ,  (?#CLSTAT 
BIT  #10000 , @#CLSTAT 
BEQ  RETURN 
MOV  #  2  ,  ERCLK 


; UNABLE  ST2  GO  BIT 
;STOP  CLOCK 

iCHECK  CLOCK  ERROR  BIT 
; NO  ERROR, GO  TO  RETURN 
;  ERROR , SET  ERCLK  2 


MOV  # 1 , DONE  ; SET  DONE  1, FINISH  SAMPLING 

RTS  PC 
.END 


oooo  o  o  ooooooooooooooo  n 


(C) 
*  *  *  * 


*  *  *  J t  * 


1 

10 


Z 

c 

c 

20 

Z 

C 

30 

Z 

C 

C 

40 

C 

C 


130 

120 

C 


79 

70 

C 

C 


Sampling  program  for  oscillating  flow, 

a**************************************,********* 

PROGRAM  DATCL 

THIS  IS  A  FORTRAN  DRIVER  FOR  COLLECTING  DATA  AND 
VRITTING  THEM  DIRECTLY  TO  A  DISK .  THIS  PROGRAM 
STARTS  TO  TAKING  DATA  WHEN  AN  OUTER  TRIGGER  SIGNAL 
OCCURS, THEN  SWEEPS  ALL  CHANNELS  AT  EACH  SAMPLING 
POINT  UNTIL  THE  REQUIRED  PERIODS  ARE  REACHED. 

DATA  ACQUISITION  RATE  IS  2000  POINTS/SEC  FOR  MULTIPLE 
CHANNELS . 

NCHS  NUMBER  OF  CHANNELS 

NPTS  NUMBER  OF  POINTS  PER  CHANNEL  PER  PERIOD 

NPRDS  NUMBER  OF  PERIODS 

DIMENSION  OF  IX  NCHS*NPTS 

BEFORE  RUN  THIS  PROGRAM,  LINK  IT  WITH  DATCLT. MAC 

********************************************************* 
PROGRAM  DATCL 
INTEGER  I  X ( 1 9  2 ) 

LOGICAL* 1  DATFI L( 1 2 ) , ERR 

VR ITE ( 7 , 10) 

FORMAT( IX, ‘TYPE  IN  THE  NAME  OF  THE  OUTPUT  FILE') 

CALL  CETSTR ( 5 , DATF I L , 1 1 , ERR ) 

IF (ERR)  GO  TO  1 


OPEN  THE  DATA  FILE, USE  DEVICE  3  TO  WRITE 

OPEN( UNIT-3 , NAME=DATF I L , TYPE= ' NEW '  , ERR=  1 
, FORM= ' UNFORMATTED '  .RECORDS  I  2  E  =  5  0 ) 


WRITE (7, 20) 

FORMAT< IX,* TYPE  IN  THE  NUMBER  OF  CHANNELS 

AND  THE  NUMBER  OF  POINTS  PER  CHANNEL  AND  NUMBER  OF  PERIODS') 
READ (5,*)  NCHS, NPTS, NPRDS 

WR ITE  <  7 , 30 ) 

FORMAT( IX, 'TYPE  IN  THE  MULTIPLE  OF  10000/SEC 
FOR  THE  SAMPLING  RATE') 

READ (5,*)  MULT 

WRITE (7, 40) 

FORMATdX,  'DATA  AQUISITION  IN  PROGRESS') 

CALL  DATCLT ( NPTS , NCHS , MULT, 32767) 

THIS  SECTION  GETS  DATA  AND  WRITES  IT  TO  DISK 
DO  120  X=l, NPRDS 
DO  130  1=1 ,NPTS*NCHS 
CALL  GETW< I X ( I ) ) 

WRITE! 3)  IX 

CALL  CLEAN 
CLOSE ( UNIT=3 ) 

WRITE ( 7 , 70  > 

FORMATdX,  'DATA  COLLECTION  COMPLETE') 


STOP 

END 
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;ft*«****A*ft**X*ft****X***********X***ft*Xftft**ft*ft*****ft***ft* 


SUBROUTINE  DATCLT . MAC 

THIS  SUBROUTINE  IS  CALLED  FROM  DATCL . FOR . 

CALL  DATCLT ( NPTS , NCHS , MULT ,327  67) 
*xxx**xxxxxx*x*x*xxxxxxx**xxxx*xxxxx****xxxxxxxxx*xx*x*x 

. SBTTL  DATCLT. MAC 
GLOBL  DATCLT 
.  GLOBL  COLOOP 
GLOBL  GETV 
.GLOBL  CLEAN 


ADSTATs 177000 
ADBUF  =ADSTAT+2 
CLSTAT= 170420 
CLBUF=CLSTAT+2 
STVEC=4  4  4 
CLVEC  =  4  4  0 
PR  7  *3  40 


THE  ADDRESS  OF  THE  A/D  BOARD 
THE  ADDRESS  OF  THE  DATA  REGISTER 
THE  ADDRESS  OF  THE  CLOCK 

ST2  INTERRUPT  VECTOR 

THE  CLOCK  INTERRUPT  VECTOR 

INTERRUPT  PRIORITY  IS  7 


P5ECT  $DATA  RV.O, 
BUFS I Z^327  67 

,  LCL  .■  RE  L  ,  CON 

NCHS: 

.WORD  0 

NPTS  : 

WORD  0 

MULT. 

.WORD  0 

LIMIT: : 

.WORD  0 

TEMP  : 

BLKV  1 

INPTR: • 

.WORD  BUF 

OUTPTR ' 

WORD  BUF 

COUNT*  : 

.WORD  0 

BUF:  : 

BLXV  BUFSIZ 

DATCLT: 

MOV  ®<R5)+,R0 

; IGNORE  NUMBER 

OF  ARGUMENTS 

MOV  ®<R5)+,N?TS 

;  GET 

NUMBER  OF 

POINTS 

MOV  @ (R5 ) +  r NCHS 

*,  GET 

NUMBER  OF 

CHANNELS 

MOV  ® ( R5 ) +  , MULT 
MOV  0<R5>+, LIMIT 

;  GET 

MULTIPLE  F 

OR  CLOCK 

MOV  tCOLOOP ,®#STVEC 
MOV  #PR7 , ®#  STVEC  f 2 
MOV  tCOLOOP ,  @#CLVEC 
MOV  #PR7 , ®#CLVEC+2 


LOAD  5T2  INTERRUPT  LOCATION 
AND  PRIORITY 
LOAD  INTERUPT  LOCATION 
AND  PRIORITY 


1$  : 


CLR  etADSTAT 
MOV  #0 ,R1 
SOB  R1  , 1  $ 


CLR  3#CLSTAT 

MOV  #BUF , INPTR 

MOV  #BUF,OUTPTR 

NEG  MULT 

MOV  MULT , @#CLBUF 

MOV  NPTS , TEMP 

MOV  #60 132 , ®#CLSTAT 

TST  SiADBUF 


; INITIALIZE  CLOCK 


} CHANGE  MULT  TO  2  * S  COMPLEMENT 
; MOV  MULT  INTO  CLOCK  COUNTER 

; ST2  START  CLOCK  t 0 0 0 0 / SEC , REPEATED  MODE 
; INTERUPT  ENABLED 


RTS  PC 


.RETURN  TO  MAIN  PROGRAM  TO 
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COLOOP: 


1  $  : 


2$ 


MOV  #O,0#ADSTAT 
BIC  #200/ 0#CLSTAT 
MOV  R0 / - ( SP ) 

MOV  R4,-<SP) 

MOV  NCHS/R4 
TSTB  0#ADSTAT 
BPL  1$ 

MOV  0#ADBUF,RO 
DEC  R4 
BE Q  3  $ 

BID  #40000, 0#ADSTAT 
JSR  PC/PUTV 
TSTB  0#ADSTAT 
BPL  25 

MOV  0#ADBUF,RO 
INCB  0# ADSTAT 
JSR  PC/PUTV 


; VA IT  FOR  INTERUPT 

; LOAD  CHANNEL  0  AND  START  CONVERSION 

; SAVE  R0 
; SAVE  R4 

;  SAVE  NUMBERS  OF  CHANNELS 
; VA I T  UNTIL  COMPLETE 

/SAVE  DATA  POINT 


; SET  SEQUENTIAL  MODE 
; SAVE  DATA  IN  BUFFER 
iTEST  IF  DONE 

/SAVE  NEXT  DATA  POINT 
; LOAD  NEXT  CHANNEL 
iSAVE  DATA  IN  BUFFER 


3$ 


EXIT; 


DEC  R4 
BNE  2$ 

DEC  TEMP 
BEQ  EXIT 
MOV  (SD  +  /R4 
MOV  <  SP )  +  , R0 
RTI 

MOV  <  SP )  +  , R4 

MOV  ( SP ) + , R0 

MOV  NPTS/TEMP 

MOV  #6 0 1 3 2 / 0#CL5TAT 

RTI 


;  DECREMENT  NCHS 

,ANY  MORE  CHANNELS 
; NO  DECREMENT  POINTS 
/EXIT  IF  DONE 


; RETURN  FROM  INTERRUPT 


; NO, BACK  FOR  ANOTHER  ST2  EVENT; 


PSECT  PUTV 

PUTV:  *  CMP  COUNT/LIMIT 

BGE  2$ 

MOV  ROf0INPTR 
INC  COUNT 
ADD  #2 , INPTR 

CMP  INPTR, #BUF+<BUFSIZ*2> 
3LO  1$ 

MOV  #BUF, INPTR 
15;  RTS  PC 

25 :  BIC  #100 ,0#CLSTAT 

RTS  PC 


GETV: : 


15  : 

CLEAN: ; 


THIS  SECTION  IS  CALLED  FROM  FORTRAN  PROGRAM  AND  MOVES  DATA  FROM 
DATAEPoInTSHTORdI^KN  CALL  STATEMENT*  TIiE  FORTRAN  PROGRAM  THEN  VRITES 

.PSECT  GETV 
CMP  #1000, COUNT 
BGE  GETV 

MOV  0OUTPTR , 02  <  R5 ) 

DEC  COUNT 
ADD  #  2 , OUTPTR 

CMP  OUTPTR ,#BUF+<BUFSIZE*2> 

BLO  1$ 

MOV  #BUF, OUTPTR 
RTS  PC 
CLR  COUNT 
CLR  0tCLSTAT 
CLRB  0#ADSTAT 
RTS  PC 
.  END 
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APPENDIX  IV.  Pulsating  Pump  Curve 


A  schematic  drawing  of  a  pulsating  pum  is  shown  in  Figure  IV. 1. 
Point  0  represents  the  shaft  of  the  motor.  Point  C,  the  position  of 
the  piston,  is  moving  along  axis  oc.  Term  Lg  is  the  stroke  length  and 
L  is  the  length  of  the  arm  which  connects  the  piston  in  one  end  and 
the  storke  in  the  another  end. 

At  any  moment,  the  position  of  the  piston  C  is  expressed  as 


x  =  L  cos  oat  + 
s 


T  2  •  2  . 

L  sm  oat 
s 


(IV. 1) 


where  oa  is  angular  frequency.  The  velocity  of  the  piston  is  then 


u 


dx 

dt 


=  L  oa  sin  oat 


-  1 


(IV. 2) 


Therefore,  the  pressure  gradient  can  be  expressed  as  -  —  4^-  «  —  and 

p  dx  dt 

plotted  in  Figure  IV. 2  for  L  =7.62  cm,  L  =  55.25  cm  and  f  =  0.325  Hz. 

s  a 

The  least  square  fit  of  a  cosine  curve  is  also  plotted  in  dashed  line.  It 

L 

s 

is  seen  that  when  is  not  negligible  the  movement  of  the  piston  starts  to 

a 


deviate  from  a  sinusoidal  variation. 
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Figure  IV. 1  Schematic  of  pulsating  pump 


Figure  IV. 2  Pressure  gradient  calculated  from  the 
movement  of  the  piston 

(L  =  7.62  cm,  L  =55.25  cm  and  f  =  0.325  Hz) 
03. 
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Symbol 

a 


a 


0 


C 

Ct 

D 

D 

m 

f 

F 

H 

I 

k 

ki 


^15k2 

K 

i 

L 


NOMENCLATURE 


Unit 

Amplitude  of  pressure  gradient  g/cm~sec 

Dimensionless  amplitude  of  flow  oscillation 

Amplitude  correction  factor  for  frequency  response 
of  the  probe;  Van  Driest  damping  constant 

Area  of  electrode 

Factor  defined  in  (3.43) 

Dimensionless  time  dependent  concentration 
defined  in  (4.10) 

Concentration 

Bulk  concentration 

Pipe  diameter 

Damping  function  in  Van  Driest  formula 
Frequency 

Faraday Ts  constant,  96,500 

Ratio  of  two  consequent  grid  sizes 

Current  flowing  in  the  electrochemical  cell 

Von  Karman  constant  0.4; 

Time  dependent  component  of  mass  transfer  coefficient  cm/sec 

Imaginary  part  of  k 

Relaxation  constant  defined  in  (3.47) 


Real  part  of  k 

Pseudosteady  approximation  of  k  cm/sec 

Coefficients  defined  in  (3.46) 

Mass  transfer  coefficient  cm/sec 

Mixing  length  cm 

Length  of  electode  (in  x  direction)  cm 


cm 


mole/cm 
mole/ cm^ 
cm 

1/ sec 

abs-coul 

q-equivalent 

amp 
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Symbol 


L 

s 


n 

e 


N 

P 

r 

ro 

Re 

s 


s 

X 


S 


Sc 

S 

X 

t 

T 


u 

* 

u 


u 


b 


u 

c 

U 


U  1 


v 

w 

W 

X 

X 


Unit 


Stroke  length  of  pulse  generator 

Number  of  electrons  transferred  in  electrochemical 
reaction 

Number  of  total  samples; 

Rate  of  mass  transfer  per  unit  area 

Pressure 

Coordinate  in  radial  direction 
Pipe  radius 
Reynolds  number 

Time  dependent  component  of  velocity 
gradient  at  the  wall 

Streamwise  velocity  gradient  fluctuation 
at  the  wall 

Velocity  gradient  at  the  wall 
Schmidt  number 

Streamwise  velocity  gradient  at  the  wall 
Time 

Dimensionless  time  defined  in  (3.53) 

Streamwise  velocity 
Friction  velocity 
Bulk  velocity  in  pipe 
Central  velocity 

Dimensionless  oscillation  velocity 

Dimensionless  time  averaged  velocity 

Radial  velocity 

Spanwise  velocity 

Width  of  the  electrode 

Coordinate  in  the  flow  direction 

Dimensionless  x  (X  =  x/L) 


cm 


mole/ cm^-sec 

,  2 
g/cm-sec 

cm 

cm 


cm/ sec 

cm/ sec 
cm/ sec 

cm/ sec 
sec 

cm/ sec 
cm/sec 
cm/ sec 
cm/ sec 


cm/ sec 
cm/sec 
cm 

cm 
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Symbol 


y 

Y 


Greek 

Symbol 


a 

0 

6 


n 

0 

A 

y 

V 

vt 

p 

X 


T(t> 

T 

W 

x  (0) 

^0 


Coordinate  normal  to  the  wall 

+  * 

Dimensionless  y.  (y  =  yu  /  v) 

*  +  1/3  +1/3 

Dimensionless  y.  (y  =  y  Sc  /L  ) 
Dimensionless  Y.  (Y  =  y/r^) 


Dif fusivity ; 

Wave  number  of  wavy  surface 

+  * 

Dimensionless  wave  number,  (a  =  av/u  ) 
Coefficient  defined  in  (4.33) 

Concentration  boundary  thickness 
Stokes  layer  thickness 
Viscous  sublayer  thickness 
Dimensionless  length  defined  in  (4.6) 

Phase  degree 

Wave  length  of  wavy  surface 
Molecular  viscosity 
Kinematic  viscosity 
Turbulent  eddy  viscosity 
Density  of  fluid 
Shear  stress 
Turbulent  shear  stress 
Wall  shear  stress 

Dimensionless  amplitude  of  wall  shear  stress, 
defined  in  (4.25) 

Coefficient  defined  in  (3.45) 

Coefficient  defined  in  (3.51) 


Unit 


cm 


2, 

cm  /sec 
1/ cm 


cm 


cm 


cm 


degree 


cm 


g/ cm-sec 

2. 

cm  /sec 
2  , 

cm  / sec 

/  3 
g/cm 

,  2 

g/cm-sec 

,  2 

g/ cm-sec 

,  2 

g/cm-sec 


Greek 

Symbol 


Unit 


<0 

Angular  frequency 

1/sec 

+ 

a o 

+  *2 

Dimensionless  to  ((a)  =  u>v/u  ) 

* 

0) 

*  +  1/3 

Dimensionless  to  (to  =  to  Sc  ) 

Q 

Frequency  parameter  (ft  =  /  rn/v  r^) 

Overline  and  brackets: 

Peoretical  component 
—  Time  average 

<  >  Phase  average 

|  |  Amplitude 

Superscripts : 

+  Nondimensionlized  with  firction  velocity  u 

and  ftinematic  viscosity  v 

(t)  Turbulent  component 
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